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Abstract 

A  25  square  kilometer  urban  area  located  in  University  Park,  Texas  was 
modeled  using  geographic  information  systems  (GIS)  software.  A  smaller 
6.25  square  kilometer  area  around  the  Southern  Methodist  University  was 
extracted  from  the  larger  area  and  analyzed  using  an  acoustic  finite- 
difference  time-domain  (FDTD)  code.  The  procedures  to  model  the  area 
using  GIS  software,  extract  required  data,  produce  the  data  files  necessary 
to  perform  the  acoustic  analyses,  and  view  the  results  are  described  in  this 
report.  This  report  is  focused  on  the  details  of  using  the  particular  software 
packages  used  in  this  study,  but  the  concepts  are  general  in  nature  and  can 
be  applied  to  other  applications  if  needed. 
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1  Overview  of  Modeling  an  Urban  Area 

1.1  Introduction 

Infrasound  signals  are  acoustic  signals  below  20Hz  and  can  be  monitored 
at  distances  of  tens  to  thousands  of  kilometers  (Hedlin  et  al.  2012).  These 
signals  are  generated  by  natural  sources  (surf,  earthquakes,  volcanoes, 
etc.),  explosions  (nuclear  and  non-nuclear),  infrastructure  sources 
(bridges,  dams,  large  buildings,  etc.),  and  human  sources  (trains,  highway 
traffic,  etc.)  (Campus  and  Christie  2010;  McKenna  et  al.  2009;  Guzas 
and  Tricys  2010).  Recent  observations  show  that  human  activity  and 
infrastructure  emit  surprisingly  ubiquitous  and  high  levels  of  infrasound. 
Because  the  wavelengths  of  infrasound  are  on  the  scale  of  tens  to  hundreds 
of  meters,  infrasonic  energy  can  propagate  tens  of  kilometers  through 
areas  of  human  activity  without  loss-of-signal  character  because  the 
wavelengths  are  transparent  to  the  audible  acoustic  obstructions,  such  as 
buildings.  Consequently,  infrasound  deployments  would  result  in  new 
human  environment-monitoring  capabilities. 

The  production  mechanisms  and  extent  to  which  different  sources  can  be 
differentiated  and  monitored  is  still  poorly  understood.  If  infrasound  is  to 
be  used  for  monitoring  the  “health”  of  a  human  environment  and  for 
identifying  anomalous  activities,  the  following  key  scientific  issues  must  be 
examined:  1)  What  are  the  typical  infrastructural  sources  of  infrasound 
and  their  levels?  2)  How  saturated  is  the  environment  of  interest  with 
infrasonic  signals  (i.e.,  do  many  signals  propagate  over  long  distances  to 
reach  a  given  sensor,  or  can  individual  sources  be  well  differentiated)?  3) 
Does  infrasound  provide  new  information  to  characterize  rapidly  evolving 
physical,  cultural,  economic,  and  military  actions  of  interest? 

Numerical  modeling  method  is  used  to  explore  signal  generation  and 
propagation  in  the  complex  urban  environment,  while  controlling  a  single 
infrasound  source  in  the  environment.  Parameterization  studies  for 
understanding  the  effect  of  the  urban  terrain  on  infrasound  signal 
propagation  and  the  effects  of  the  short  propagation  distances  related  to 
local  metrological  effects  can  be  completed  while  controlling  the  sources 
and  other  variables  in  the  environment.  As  the  parameterization  studies 
are  completed,  base  cases  are  validated  with  real-world  data. 
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1.2  Purpose  of  report 

This  report  describes  the  modeling  of  the  structures  and  landscape  in  the 
urban  area  with  sufficient  resolution  and  accuracy  to  allow  the  subsequent 
analysis  with  an  acoustic  finite-difference  time-domain  (FDTD)  code.  The 
purpose  of  the  report  is  to  lead  the  reader  through  the  steps  involved  in 
gathering  the  required  geographical  data,  modeling  the  3-dimensional  (3-D) 
urban  area  of  interest  in  geographical  information  systems  (GIS)  software, 
production  of  needed  data  files,  analysis  requirements,  and  viewing  of 
results.  This  report  details  techniques  and  methods  to  evaluate  the  princi¬ 
ples  of  urban  infrasound  propagation.  This  report  is  focused  on  the  details 
of  using  various  software  packages  to  perform  the  modeling,  but  some  of  the 
concepts  are  general  in  nature  and  can  be  applied  to  other  applications  if 
needed.  The  software  packages  used  are  discussed  in  Section  1.4. 

1.3  Description  of  urban  area  modeied 

A  25  km2  urban  area  was  modeled  in  support  of  infrasound  studies  related 
to  acoustic  propagation  through  urban  environments.  The  area  was  in 
University  Park,  Texas,  around  the  Southern  Methodist  University,  as 
shown  by  the  red  outlined  area  in  Figure  1.1.  The  area  measured  5  km  on 
each  side.  A  smaller  6.25  km^  area  was  used  for  analysis  purposes,  as  shown 
by  the  green  outlined  area  in  Figure  1.2.  This  area  was  2.5  km  on  each  side. 
The  location  of  the  lower-left  corner  of  the  areas  shown  in  Figures  1.1  and 
1.2  are  given  in  the  Universal  Transverse  Mercator  (UTM)  coordinate 
system. 

The  work  required  the  synthesis  of  GIS  data  from  various  sources  to 
produce  a  3-D  model  of  the  urban  area.  The  main  GIS  data  used  pertained 
to  defining  the  topography,  the  paved  areas  (e.g.  road,  highways,  and 
parking  lots),  and  the  buildings.  The  orange  areas  shown  in  Figures  i.i  and 
1.2  are  buildings.  The  gray  areas  shown  in  Figures  1.1  and  1.2  are  the  paved 
areas.  The  GIS  data  is  more  fully  described  in  Chapter  2. 

Analyses  were  performed  on  the  area  shown  in  Figures  1.2  and  1.3  using  an 
acoustic  FDTD  code  running  on  the  ERDC  Major  Shared  Resource  Center 
supercomputers . 
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(706301.196,  3634845.629) 


Figure  1.1.  25  km^  area  in  University  Park,  Texas  with  smaller 
(green)  analysis  area. 


(704837.812,  3633790.296) 

Figure  1.2.  2.5  km^area  around  Southern  Methodist  University. 
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Figure  1.3.  3-D  view  of  area  around  Southern  Methodist  University. 


1.4  Discussion  of  programs  used  and  steps  to  produce  urban  data 
sets  and  perform  anaiyses 

Several  commercial  and  custom  programs  were  used  to  produce  the  urban 
data  files  and  to  conduct  acoustic  analyses  of  the  6.25  km^area  shown  in 
Figure  1.2.  A  short  discussion  of  each  program  follows  to  explain  how  the 
programs  were  used  to  process  the  GIS  urban  data  sets.  The  GIS  data  sets 
will  be  explained  more  fully  in  Chapter  2.  The  process  steps  are  shown  in 
Figure  1.4  below. 

Two  GIS  programs  were  used  to  manipulate  the  data  sets  used  to  model  the 
urban  area.  Global  Mapper  (version  14.0)  from  Blue  Marble  Geographies 
(2013)  and  ArcMap  (version  10.0)  from  the  ArcGIS  product  from 
Environmental  Systems  Research  Institute  (ESRI 2013). 

Global  Mapper  was  used  to: 

•  Convert  all  data  to  a  common  projected  coordinate  system. 

•  Trim  large  data  sets  to  the  extents  of  the  desired  analysis  area. 

•  Modify  the  elevation  of  the  topography  under  bridges. 

•  Extract  certain  features  from  larger  data  sets  based  on  location,  size,  or 
height. 

•  View  2-D  buildings  in  3-D  given  a  height  attribute. 

•  Output  ASCII  vector  files  defining  the  needed  data  sets  for  additional 
processing  by  other  software. 
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Figure  1.4.  Steps  in  modeling  and  running  an  acoustic  analysis  of  an  urban  area. 
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The  output  of  ASCII  files  by  Global  Mapper  made  for  an  easy  path  to  read 
the  data  and  perform  other  custom  operations. 

ArcMap  was  used  to: 

•  Perform  custom  queries  of  data  to  segment  the  data  by  desired  ranges 
of  attributes. 

•  Perform  queries  based  on  containment  of  one  feature  in  another. 

•  Make  needed  data  sets  by  performing  relational  database  joins  on  the 
information  contained  in  the  geodatabase. 
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•  Aggregate  buildings  into  larger  polygons  that  represented  buildings 
with  a  common  range  of  a  desired  attribute  (i.e.,  height).  This  also 
effectively  segmented  the  buildings  based  on  function,  such  as 
residential  and  commercial. 

•  Compute  weighted  averages  of  heights  of  buildings  contained  with 
larger  polygons 

•  Make  data  sets  based  on  the  proximity  of  a  feature  set  (e.g.,  buildings) 
to  other  geographic  features  (e.g.,  roads). 

•  Output  features  in  the  ESRI  shapefile  format  for  reading  into  Global 
Mapper. 

ArcMap  was  used  to  aggregate  the  buildings  in  the  urban  environment 
into  larger  areas  based  on  the  proximity  of  buildings  to  one  another  and  in 
computing  weighted  averages  of  heights  of  buildings  within  desired 
polygons.  This  is  further  explained  in  Chapter  3. 

A  piece  of  custom  software  termed  the  Urban  Modeler  was  used  to  read 
the  sets  of  data  in  the  ASCII  format  output  by  Global  Mapper.  Urban 
Modeler  was  then  used  to: 

•  Convert  the  units  of  vertical  measurement  of  attributes  to  a  common 
unit  of  meters. 

•  Determine  if  islands  exist  in  the  data  and  write  the  island  information 
out  separately  from  the  feature  information. 

•  Break  large  sets  of  data  (e.g.,  thousands  of  buildings)  into  many 
smaller  sets  of  data  to  improve  processing  later  in  MATLAB. 

•  Produce  data  files  for  buildings  and  area  features  (i.e.,  paved  features 
consisting  of  streets,  highways,  parking  lots,  driveways,  and  medians). 

•  Write  GIS  information  out  to  ASCII  comma  separated  variable  (CSV) 
files  for  use  with  Microsoft  Excel. 

Islands  represent  open  spaces  in  a  larger  area.  That  is,  a  building  may  exist 
with  an  interior  courtyard,  which  can  be  modeled  using  an  island.  Paved 
areas,  such  as  highways,  may  have  interior  open  spaces  that  represent  the 
medians.  The  network  of  streets  may  be  represented  as  large  solid  areas 
with  islands  representing  the  city  blocks.  The  GIS  programs  will  display 
the  features  and  the  islands  will  cause  the  large  solid  areas  to  appear  to 
have  open  spaces  in  them.  The  islands  must  be  accounted  for  in  the 
modeling  process.  This  is  discussed  further  in  Chapter  4. 
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The  CSV  files  produced  by  the  Urban  Modeler  program  define  the 
structures  (buildings  and  paved  areas)  that  sit  on  the  topography.  These 
files  were  read  into  Microsoft  Excel  (version  2007)  and  then  inserted  into 
other  Excel  spreadsheets  that  contained  the  definition  of  the  problem 
domain.  This  information  consisted  of  the  definition  of  the  analysis  grid 
size,  default  material  property  for  the  ground,  air  layers  and  properties, 
topography,  and  the  structures  that  reside  on  the  topography.  The  use  of 
the  Excel  spreadsheets  is  explained  more  in  Chapter  5. 

The  MATLAB  (version  2012b)  program  from  MathWorks  (2013)  was  used 
to  run  custom  scripts  that  processed  the  information  contained  in  the 
Excel  spreadsheets  and  produced  output  used  to  compile  and  run  the  3-D 
acoustic  FDTD  program  called  PSTOP3D.  The  use  of  the  MATLAB  scripts 
is  explained  in  Chapter  6. 

PSTOP3D  (version  AS217),  which  roughly  stands  for  parallel  stair-casing 
topography  in  3-D,  was  used  to  perform  an  acoustical  analysis  of  the  3-D 
data  representing  the  urban  environment.  PSTOP3D  is  a  FDTD  code 
produced  by  the  ERDC  Cold  Regions  Research  Laboratory  (Ketcham 
2006;  Ketcham  et  al.  2005,  2007,  2008;  Parker  et  al.  2007).  The 
procedures  to  compile  and  run  PSTOP3D  are  explained  in  Chapter  7. 

MATLAB  was  also  used  to  post-process  the  results  to  produce  various 
plots  and  animations.  The  post-processing  is  discussed  in  Chapter  8. 
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2  GIS  Data  Sets  Used  to  Model  the  Urban 
Area 

2.1  Sources  for  GIS  data 

The  information  contained  in  this  chapter  corresponds  with  Step  i  in 
Figure  1.4.  There  are  many  commercial  and  government  sources  of  GIS 
data.  Some  of  the  government  sources  that  were  investigated  in  this  study 
were  the: 

•  Army  Geospatial  Center  (AGC) 

•  National  Geospatial  Intelligence  Agency  (NGA) 

•  United  States  Geological  Service  (USGS) 

•  United  States  Departments  of  Agriculture  (USDA) 

•  United  States  Census  Bureau 

These  agencies  have  information  consisting  of: 

•  Digital  elevation  data  at  varying  resolutions 

•  Photographic  data  consisting  of  satellite  and  aerial  imagery 

•  Land  cover  and  land  use  data 

•  Road  and  highway  data 

The  two  state  agencies  that  were  used  to  collect  data  were  the  North 
Central  Texas  Council  of  Governments  (NCTCOG)  and  the  Texas  Natural 
Resources  Information  System  (TNRIS). 

The  data  sets  of  interest  defined  the: 

•  Topography 

•  Individual  building  footprints  with  a  height  attribute 

•  Areal  extents  of  paved  areas  such  as  highways,  roads,  streets,  parking 
lots,  and  medians 

•  Areal  extents  of  water  features 

•  Areal  extents  of  trees  or  forests  with  height  information 
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2.2  Topography  data 

Sources  given  in  Section  2.1  for  topographical  data  were  investigated  and 
several  sets  of  digital  elevation  data  were  collected.  The  data  set  determined 
to  best  represent  the  topography  of  the  given  area  was  obtained  from  TNRIS 
(2013). 

The  topography  was  modeled  by  using  light  detection  and  ranging 
(LIDAR)  data  purchased  from  TNRIS  which  is  part  of  the  Texas  Water 
Development  Board.  The  LIDAR  data  was  flown  during  2009  and 
consisted  of  data  points  spaced  at  1  meter. 

LIDAR  systems  are  typically  mounted  on  an  aircraft  and  use  a  laser  to 
densely  sample  the  surface  of  the  earth  and  produce  highly  accurate 
terrain  models. 

The  data  set  was  provided  in  several  formats  including  LIDAR  laser  (LAS) 
files  and  files  in  the  United  States  Geological  Service  (USGS)  digital 
elevation  model  (DEM)  format.  The  LAS  files  contain  information  on  the 
classification  of  the  LIDAR  points,  the  return  value,  and  the  intensity 
value.  Points  are  classified  as  objects  such  as  trees,  water,  and  buildings. 
The  last  return  value  represents  the  bare  earth  model. 

The  DEM  files  contained  the  bare  earth  model  that  could  also  be  extracted 
from  the  LAS  files.  This  study  used  the  DEM  files  to  input  the  topography 
into  Global  Mapper.  The  topography  defined  by  the  LIDAR  data  is  shown 
in  Figure  2.1. 

2.3  Aerial  orthophotography 

2009  6  in.  aerial  orthophotography  was  obtained  from  NCTCOG  in  the 
GeoTIFF  image  format.  The  imagery  was  mainly  used  to  provide  a  nice 
overlay  for  the  topography  when  producing  displays  of  the  3-D  modeled 
area. 


Additional  orthophotography  was  obtained  from  the  USDA  National 
Agriculture  Imagery  Program  (NAIP  2013)  obtained  in  2012.  This  imagery 
was  used  to  add  a  custom  building  to  the  building  data  set. 


ERDC  TR-15-5 


10 


Figure  2.1.  Topography  of  6.25  km^area. 


2.4  Building  data 

Data  were  sought  that  would  provide  a  definition  of  the  building  footprints 
with  an  associated  height  attribute.  This  would  be  enough  information  to 
construct  an  extruded  building  with  the  appropriate  footprint  and  a 
representative  height  to  adequately  represent  the  building  in  an  analysis. 

Several  sources  of  data  were  examined  as  described  in  Section  2.1.  The 
AGC  had  information  from  NGA  from  a  project  called  the  133  Cities 
Project.  This  project  collected  LIDAR  data  on  urban  areas  and  processed 
the  data  to  produce  a  bare  earth  model,  footprints  of  buildings  with  a 
height  attribute,  and  definition  of  trees  and  forested  areas.  These  data 
would  have  been  perfect  except  that  data  for  only  a  portion  of  the  desired 
area  were  available. 

Planimetric  data  sets  from  the  NCTCOG  were  available  for  the  entire  area. 
Planimetric  data  consist  of  2-D  representations  of  objects  as  seen  from 
aerial  photography.  The  features  include  roads,  building  footprints, 
sidewalks,  trails,  rivers,  lakes,  etc.  These  features  are  often  digitized  from 
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orthorectified  aerial  photography  producing  images  that  have  been 
adjusted  for  topographic  relief,  lens  distortion,  and  camera  tilt.  The  images 
may  be  used  to  accurately  measure  distances.  Unfortunately,  the  building 
footprints  contained  in  the  planimetric  data  did  not  have  height  attributes. 

Eventually,  a  contract  was  undertaken  with  CyberCity  3D  (2012),  a 
commercial  company  located  in  El  Segundo,  CA,  to  construct  the 
buildings.  The  company  specializes  in  the  production  of  accurate  3-D 
models  of  cities  including  3-D  buildings,  streets,  trees,  and  urban  objects 
that  can  be  used  for  various  visualization  of  planning  functions. 

CyberCity  3D  used  6  in.  stereoscopic  aerial  photography  to  construct  3-D 
models  of  the  buildings.  These  models  were  quite  detailed  and  included  the 
shape  of  the  roof  (e.g.  flat,  hip,  gable,  gambrel)  and  additional  shapes,  such 
as  dormers.  CyberCity  3D  processed  the  data  using  their  in-house  custom 
software. 

A  set  of  2-D  buildings  were  also  produced  that  included  the  footprints  of 
the  buildings  and  a  height  attribute.  The  height  of  the  building  was  defined 
as  the  lower  edge  of  the  building  up  to  the  mid-height  of  the  roof.  For  a  flat 
roof,  the  measurement  is  from  the  lowest  bottom  edge  of  the  building  to 
the  top  of  the  roof.  For  a  roof  that  is  not  flat,  such  as  a  hip  or  gable  roof, 
the  measurement  is  from  the  lowest  bottom  edge  of  the  building  up  to  the 
mid-height  of  the  roof.  The  height  measurement  was  provided  in  feet. 

The  data  set  was  provided  as  ESRI  shapefiles  and  AutoCAD  DXF  files.  The 
shapefiles  were  used  for  the  2-D  footprints  while  the  AutoCAD  DXF  files 
were  used  to  represent  the  3-D  features. 

The  2-D  definitions  of  the  buildings  were  used  in  this  study  to  determine 
the  effect  of  buildings  on  infrasound  propagation.  The  footprints  of  the 
buildings  are  shown  in  Figures  1.1  and  1.2  as  orange  areas.  The  larger  area 
shown  in  Figure  1.1  consisted  of  54,225  footprints.  The  smaller  area  shown 
in  Figure  1.2  consisted  of  12,568  footprints.  This  is  not  the  total  number  of 
buildings  because  a  single  building  may  have  multiple  footprints  defining 
different  levels  or  pieces  of  geometry  such  as  dormers.  The  actual  number 
of  buildings  used  in  the  analysis  for  the  area  shown  in  Figure  1.2  is  11,059. 
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2.5  Paved  areas  data 

Data  were  needed  to  define  paved  areas  such  as  roads,  streets,  highways, 
parking  lots,  and  medians.  Data  defining  the  roads,  streets,  and  highways 
were  found  from  the  U.S.  Census  Bureau  as  TIGER/Line  (Topologically 
Integrated  Geographic  Encoding  and  Referencing)  files.  These  files  were 
ESRI  shapefiles  and  defined  the  roads,  streets,  and  highways  in  sufficient 
detail  topologically.  Attributes  were  provided  that  gave  the  feature  name 
and  the  type  of  the  feature  such  as  residential  road.  The  problem  with  the 
data  was  that  the  roads  did  not  have  a  width  specified.  Therefore,  a  width 
would  have  to  be  assumed  based  upon  the  feature  type  specified. 

Planimetric  data  from  the  NCTCOG  were  available  for  the  entire  area  for 
2007.  Planimetric  data  consist  of  2-D  representations  of  objects  as  seen 
from  aerial  photography.  The  features  include  items  such  as  roads, 
building  footprints,  sidewalks,  trails,  rivers,  lakes,  bridges,  etc.,  which 
were  traced  using  digital  orthorectified  aerial  photography,  thereby 
resulting  in  the  accurate  areal  definition  of  the  objects. 

All  visible  paved  parking  lots  large  enough  to  accommodate  five 
automobiles  were  digitized.  No  gravel  or  dirt  parking  lots  were  digitized. 
The  accuracy  of  all  digitized  features  matches  the  accuracy  of  the  aerial 
photography  used.  Aerial  photography  from  NCTCOG  meets  National  Map 
Accuracy  Standards  at  a  map  scale  of  1  in.  =  200  ft  (1:2400). 

The  planimetric  data  set  was  useful  because  it  defined  the  areal  extents  of 
the  roads,  streets,  highways,  and  bridges.  The  data  set  also  provided 
information  for  parking  lots,  medians,  and  sidewalks.  Since  the  areal 
extents  were  defined,  a  width  did  not  have  to  be  assumed  for  the  roads, 
streets,  or  highways.  The  sidewalk  information  was  not  used  in  the 
analyses,  but  was  modeled  by  using  buffer  zones  of  a  specified  size  that 
corresponded  to  an  average  sidewalk  width.  This  process  is  explained  in 
Chapter  3.  The  planimetric  files  were  provided  as  ESRI  shapefiles. 

For  the  area  defined  in  Figure  1.1,  there  were  2,155  paved-area  objects.  The 
trimmed  objects  for  the  smaller  area  shown  in  Figure  1.2  consisted  of  648 
objects.  The  objects  could  cover  very  large  areal  extents  and  this  caused 
some  problems  with  the  data  processing.  To  alleviate  the  problems,  larger 
objects  were  split  into  smaller  objects  for  processing  which  resulted  in  a 
total  of  947  objects  being  processed.  This  is  explained  in  Chapter  3. 
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2.6  Trees  and  water  features  data 

Other  objects  of  interest  in  the  urban  environment  are  trees  (individual), 
forests  (conglomeration  of  individual  trees),  and  water  features  such  as 
lakes,  streams,  and  rivers. 

For  the  analyses  performed,  trees  and  water  features  were  not  included. 
The  data  sources  discussed  in  Section  2.1  did  contain  some  information 
concerning  these  features.  For  trees  and  forests,  the  data  sets  from  AGO 
were  defined  using  LIDAR  data  and  included  a  height  attribute.  The  area 
under  consideration  did  not  have  complete  coverage.  The  LIDAR  data 
from  TNRIS  did  include  point  feature  definitions  of  trees  that  could  be 
used  to  model  the  tree  coverage.  This  was  not  done  in  this  study. 

The  water  features  were  not  adequately  defined  in  any  of  the  data  sets 
obtained.  Water  features  were  not  included  in  the  initial  analyses,  but  were 
modeled  by  simply  tracing  water  features  shown  on  2009  6  in. 
orthophotography  obtained  from  NCTCOG.  This  allowed  the  width  of  the 
features  to  be  accurately  modeled. 

2.7  Projected  coordinate  systems 

All  GIS  data  sets  are  defined  using  a  geographic  or  projected  coordinate 
system.  An  example  of  a  geographic  coordinate  system  is  the  typical 
latitude-longitude  that  defines  a  point  on  the  earth  using  degrees  of 
latitude  and  longitude.  A  projected  coordinate  system  maps  the  points  on 
a  sphere  to  a  flat  map.  There  are  many  projected  coordinate  systems  and 
you  must  know  what  system  is  used  to  define  your  data.  The  projected 
coordinate  system  information  may  be  supplied  with  the  GIS  information 
as  metadata  files,  as  projection  files,  or  from  information  provided  by  the 
GIS  data  provider.  For  the  data  described  in  the  report,  the  projected 
coordinate  systems  were  supplied  with  the  data  received  from  the  various 
data  sources. 

Even  if  the  projected  coordinate  system  of  your  data  is  known,  attributes 
supplied  with  the  data  may  be  in  different  units.  For  example,  the  data  set 
defining  the  areal  footprints  of  buildings  in  this  study  was  in  meters,  but 
the  attribute  defining  the  heights  of  the  buildings  was  in  feet. 

For  this  study,  the  data  sets  that  were  obtained  were  not  all  in  the  same 
projected  coordinate  system.  Therefore,  all  data  sets  were  re-projected  to 
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the  UTM  coordinate  system  which  Global  Mapper  does  on  the  fly.  That  is, 
the  first  data  set  read  into  Global  Mapper  sets  the  default  coordinate 
system.  All  other  data  sets  read  in  afterwards  are  re-projected  to  this  default 
coordinate  system. 

The  UTM  coordinate  system  uses  a  2-D  Cartesian  coordinate  system  to 
give  locations  on  the  surface  of  the  earth.  The  horizontal  and  vertical 
locations  are  given  in  meters.  The  UTM  system  divides  the  earth  into  sixty 
zones  that  are  each  composed  of  six  degrees  of  longitude.  This  system 
must  also  reference  a  vertical  datum.  The  one  used  in  this  study  was  the 
World  Geodetic  System  (WGS)  1984  datum.  All  data  in  this  study  fell 
within  UTM  zone  14. 

Projecting  data  to  the  UTM  coordinate  system  provides  a  2-D  Cartesian 
coordinate  system,  in  which  accurate  linear  measurements  can  be 
performed. 
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3  Process  GIS  Data 

3.1  Overview 

The  information  in  this  chapter  pertains  to  Steps  2  and  3  in  Figure  1.4  of 
the  process  map.  The  processing  of  the  GIS  data  was  performed  by  the 
program  Global  Mapper  from  Blue  Marble  Corporation  and  ArcMap  from 
ESRI.  These  programs  were  used  to  convert  the  data  sets  into  one 
common  geographic  coordinate  system,  trim  the  data  sets  to  the  needed 
areal  extents,  extract  certain  features  from  the  data,  and  aggregate  the 
building  data  into  larger  more  general  areas.  The  specific  steps  to 
accomplish  these  tasks  are  discussed  within  this  chapter. 

Although  most  of  the  tasks  could  be  performed  with  either  program,  some 
tasks  were  easier  to  perform  with  Global  Mapper  and  others  with  ArcMap. 
The  discussion  in  this  chapter  will  show  the  steps  required  to  perform  a 
task  using  either  Global  Mapper  or  ArcMap,  but  not  both.  There  are  some 
tasks,  such  as  the  aggregation  of  buildings,  which  could  only  be  performed 
in  one  program. 

Global  Mapper  was  mainly  used  to  perform  the  manipulations  and  display 
of  the  data  (both  2-D  and  3-D),  while  ArcMap  was  used  to  perform  a  few 
specific  tasks. 

3.2  Define  topography 

3.2.1  Read  topography  into  Global  Mapper 

Global  Mapper  was  used  to  read  the  topography  information  obtained 
from  TNRIS  in  the  USGS  DEM  format.  The  topography  for  the  6.25  km^ 
area  is  shown  in  Figure  2.1. 

Data  are  read  into  Global  Mapper  by  opening  the  specific  file  using  the 
File  menu  option  or  dragging  the  file  from  Windows  Explorer.  If  the  file 
contains  a  definition  of  the  projected  coordinate  system  used.  Global 
Mapper  will  automatically  read  and  display  the  data.  If  a  projected 
coordinate  system  is  not  found,  a  dialog  box  will  display  asking  the  user  to 
verify  the  appropriate  projected  coordinate  system  as  shown  in  Figure  3.1. 
The  user  would  select  the  appropriate  projected  coordinate  system  that  the 
data  use  and  select  the  OK  button. 
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Figure  3.1.  Selection  of  geographic  coordinate 
system. 


3.2.2  Create  clip  region  for  topography 

Topography  was  obtained  for  an  area  slightly  larger  than  the  25  km^  area 
shown  in  Figures  1.1  and  3.2.  The  red  outlined  area  in  Figure  3.2  is  the  clip 
area  for  the  topography  for  the  smaller  area  of  interest.  The  green  outlined 
area  is  the  area  of  interest  used  for  the  analysis.  A  larger  topography  is  used 
for  the  analysis  area  to  give  the  MATLAB  code  information  to  interpolate 
elevations  at  the  edges  of  the  analysis  area.  The  topography  for  the  analysis 
area  was  clipped  from  the  larger  area  using  Global  Mapper  as  follows: 

1.  Create  a  box  that  surrounds  the  area  of  interest.  In  Global  Mapper,  a  box  is 
created  by  selecting  the  Digitizer  Tool  button  and  then  selecting  the 
Add  Rectangular/Square  Area  button  as  shown  in  Figure  3.3. 

2.  Left-click  a  starting  point  and  drag  out  the  square  or  rectangular  region 
desired.  If  the  SHIFT  key  is  held  down,  a  square  is  created  when  dragging 
the  mouse. 

3.  When  the  box  is  created,  a  dialog  pops  up  to  name  the  feature  as  shown  in 
Figure  3.4.  Type  a  name  and  click  OK.  This  feature  will  be  used  when 
saving  the  topography  data  as  described  in  Section  3.2.3.  Another  feature 
can  be  created  to  clip  data  to  the  analysis  region. 
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Figure  3.2.  Topography  and  clip  areas  for  analysis  area. 


Figure  3.3.  Adding  an  area  to  provide  a  clip  boundary. 

Grobdt  hAapper  vl4.0  (b09l2L2)  [64-bit]  -  REGISTERED  (]Univvi^ityPark-LarigeArtd4isslivintnt£^^ 
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Figure  3.4.  Name  the  clip  feature. 


3.2.3  Create  GeoTIFF  for  MATLAB 

The  MATLAB  processing  routines  require  that  the  digital  elevation  model 

data  be  saved  as  a  32-bit  GeoTIFF  file  with  the  geographic  coordinate 

system  set  as  latitude-longitude.  This  GeoTIFF  file  can  be  created  in 

Global  Mapper  as  follows: 

1.  Turn  on  the  Overlay  Control  Center  by  left-clicking  the  icon  shown  in 
Figure  3.5.  Select  objects  to  display  by  left-clicking  the  object  to  display  a 
check  mark.  Select  all  the  topography  and  any  edits  to  the  topography, 
such  as,  flattened  areas  under  bridges  as  described  in  Section  3.5.2.  Turn 
on  a  clip  outline  if  desired  as  discussed  in  3.2.2.  Turn  off  all  other  objects. 
Flattening  the  terrain  under  bridges  is  described  in  Section  3.5.2. 

2.  From  the  Tools  menu  option,  select  the  Configure  option  and  change 
the  display  projection  to  Geographic  Oatitude/longitude)  as  shown  in 
Figure  3.6  and  click  OK  The  GeoTIFF  file  needs  to  be  in  geographic 
coordinates  versus  UTM  coordinates  for  use  in  the  MATLAB  routines.  The 
MATLAB  routines  automatically  apply  a  display  projection  to  convert 
geographic  coordinates  to  UTM. 

3.  Select  the  Feature  Info  Tool  as  shown  in  Figure  3.7  and  click  on  the  clip 
feature  to  select  it. 

4.  Select  File  >  Export  >  Export  Elevation  Grid  Format  as  shown  in 
Figure  3.8. 

5.  Select  the  option  Elevation  (32-bit  floating  point  samples)  as  shown 
in  Figure  3.9. 
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6.  Click  on  the  Export  Bounds  tab  shown  in  Figure  3.10  and  click  on  the 
Crop  to  Selected  Area  Feature(s)  option.  Click  OK  and  name  the 
export  file. 


Figure  3.5.  View  the  overiay  controi  center. 


Gloiul  Mappei  vl4.Q  [b09L2L2i  [64-bit)  -  REGISTERED  (UnivtfiilyPjifc-JSOOa'Mja-faf  repcttgow) 


jp' 


zl*l 


Seiect  overlay  control  center 


Figure  3.6.  Changing  projection 
system  to  geographic. 


Figure  3.7.  Select  feature  info  tool. 

1^  Global  Mapper  vl4.0  (b091212)  164-bit)  -  REGISTERED  (UniversityPark-2500area-for  reportgmw) 


& 

^  <3 

, 

^  f  4  & 

Feature  info  tool 

I 
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Figure  3.8.  Exporting  the  DEM  data  as  a  GeoTiFF. 


Figure  3.9.  GeoTiFF  options  in  Giobai 
Mapper. 


GeoTiFF  Options  |  Gridding  |  Export  Bounds  | 


-Fite  Type - 

<~  8-bit  Paiette  image 

F'  24-bit  RGB  (Fuii  Color,  May  Create  Large  Files) 

C"  Black  and  White  (1  bit  per  pixel) 

C"  Multi-Band  l|  16  |  -bits  per  Band)  [5  Bands 

Bevation  (16  bit  integi^sarTiptes) 

Sevabon  (32  bit  floabng  point  samptes)j 

Vertical  Units  |  Meters  ▼  | 

Resampling  |  Default  (Resample  if  Needed)  | 

-  Sample  Spacing/Scate - 

X^s;  1 1 .58579906235476e-006  arc  degrees 

Y-axis;  1 1 .2971 741 28351 67e-006  arc  degrees 

W  Always  Generate  Square  Pixels 

If  you  wish  to  char>ge  the  ground  units  that  the  spacing  is 
specified  in.  you  need  to  change  the  current  protection  by 
going  to  Con^-> Projection. 

Click  Here  to  Calculate  Spacing  in  Other  Units...  | 

I  Export  at  the  Fixed  Scale  1 :  [o 


r”  Save  Map  Layout  (Scale/M argins/Grid/Legend/etc.) 
I  Save  Vector  Data  if  Displayed 
I  Interpolate  to  Fill  SmaN  Gaps  in  Data 
1“  Generate  TFW  (World)  File 
V  Generate  PRJ  Rte 

1“  ADVANCED:  Don\  Wrte  GeoTiFF  Header 


1  ox  1 

1  Cancel  | 

Apply  1 

Help 
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Figure  3.10.  Select  the  export  bounds  to 
trim  data  to  a  specified  extent. 


3.3  Model  buildings 

The  building  data  set  provided  by  CyberCity  3D  was  processed  to  produce 

several  building  data  sets  used  for  analysis.  The  building  data  set  was 

processed  to: 

1.  Construct  a  building  data  set  that  eliminated  buildings  that  were  under 
1  m2  in  area  and  2  m  in  height.  This  eliminated  small  inconsequential 
buildings  or  pieces  of  buildings  and  aided  in  the  computation  of  the  finite 
difference  grid. 

2.  Further  divide  the  buildings  into  two  data  sets  that  consisted  of  buildings 
less  than  or  equal  to  51  ft  and  greater  than  51  ft  in  height. 

3.  For  the  building  data  set  that  was  less  than  51  ft  in  height,  the  data  set  was 
further  divided  into  buildings  that  were  within  200  meters  of  U.S. 
Highway  75  and  buildings  that  were  further  than  200  meters  from  U.S. 
Highway  75. 

4.  Both  building  data  sets  that  were  less  than  or  equal  to  51  ft  in  height  were 
then  aggregated  based  on  an  aggregation  distance  of  25  m.  This  process 
combined  buildings  that  were  within  the  aggregation  distance  into  larger 
buildings. 
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5.  The  height  of  the  aggregated  building  data  sets  was  computed  as  a 
weighted  average  of  the  individual  buildings  contained  in  the  larger 
aggregated  polygons. 

The  above  processing  results  in  five  data  sets  were  used  to  perform  the 

acoustic  analyses.  The  data  sets  consisted  of: 

1.  Buildings  that  were  greater  than  51  ft  in  height  (orange  buildings  in 
Figure  3.11) 

2.  Buildings  that  were  200  m  away  from  U.S.  Highway  75  and  less  than  or 
equal  to  51  ft  in  height  (blue  buildings  in  Figure  3.12). 

3.  Buildings  that  were  within  200  m  of  U.S.  Highway  75  and  less  than  or 
equal  to  51  ft  in  height  (pink  buildings  in  Figure  3.12). 

4.  Aggregated  buildings  that  were  200  m  away  from  U.S.  Highway  75  and  less 
than  or  equal  to  51  ft  in  height  (purple  buildings  in  Figure  3.13). 

5.  Aggregated  buildings  that  were  within  200  m  of  U.S.  Highway  75  and  less 
than  or  equal  to  51  ft  in  height  (green  buildings  in  Figure  3.13). 


Figure  3.11.  Buildings  greater  than  (orange)  and 
less  than  or  equal  to  (blue)  51  ft  in  height. 
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Figure  3.12.  Building  greater  than  (blue)  and  within  (pink) 
200  m  of  U.S.  Highway  75. 


Figure  3.13.  Aggregated  buildings  greater  than  (purple)  and  within 
(green)  200  m  of  U.S.  Highway  75. 
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The  following  sections  will  provide  the  processing  information  necessary 
to  produce  these  data  sets. 


3.3.1  Read  buildings  into  Global  Mapper 

Global  Mapper  was  used  to  read  2-D  building  data  produced  by  CyberCity 
3D  in  the  ESRI  shapefile  format.  The  buildings  for  the  25  km^  area  are 
shown  in  Figure  3.14  in  the  red  outlined  area  and  are  colored  orange.  The 
buildings  from  this  larger  area  may  be  trimmed  to  the  smaller  6.25  km^ 
area  that  is  outlined  in  green  as  detailed  in  the  following  section. 


Figure  3.14.  Buildings  for  large  area  (red)  and  area  of  interest  (green). 


Area  of  interest 


(704837.812,  3633790.296) 


Data  are  read  into  Global  Mapper  by  opening  the  specific  file  using  the 
File  menu  option  or  dragging  the  file  from  Windows  Explorer.  If  the  file 
contains  a  definition  of  the  projected  coordinate  system  used,  Global 
Mapper  will  automatically  read  and  display  the  data.  If  a  projected 
coordinate  system  is  not  found,  a  dialog  box  will  display  asking  the  user  to 
verify  the  appropriate  projected  coordinate  system  as  shown  in  Figure  3.1. 
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The  user  would  select  the  appropriate  projected  coordinate  system  that  the 
data  set  uses  and  select  the  OK  button. 

3.3.2  Clip  buildings  to  desired  area 

The  building  data  must  be  clipped  (or  trimmed)  to  the  area  used  for  the 
analysis.  This  clipping  operation  will  clip  features  located  on  the  edges  of  a 
specified  clip  area  in  order  to  have  edges  that  correspond  with  the  clip 
feature.  All  features  inside  the  clip  area  will  not  be  affected  and  all  features 
totally  outside  of  the  clip  area  will  be  discarded.  This  produces  a  set  of 
buildings  that  exactly  fits  inside  the  clip  area. 

The  construction  of  a  clip  area  is  explained  in  Section  3.2.2.  After  the  clip 
area  is  constructed,  the  steps  to  clip  and  save  the  data  in  Global  Mapper 
are  as  follows: 

1.  Turn  on  the  Overlay  Control  Center  by  left-clicking  the  icon  shown  in 
Figure  3.5.  Select  objects  to  display  by  left-clicking  the  object  to  display  a 
check  mark.  Turn  all  layers  off  except  for  the  building  data  that  is  to  be 
clipped.  Turn  on  the  clip  area  as  discussed  in  3.2.2. 

2.  If  the  display  projection  is  not  already  UTM,  it  must  be  changed  to  UTM.  The 
MATLAB  routines  can  read  data  in  the  UTM  coordinate  system  and  this 
makes  inputting  the  buildings  easier  than  using  a  local  coordinate  system. 
From  the  Tools  menu  option,  select  the  Configure  option  and  change  the 
display  projection  to  UTM  as  shown  in  Figure  3.15  and  click  OK 

3.  Select  the  Feature  Info  Tool  as  shown  in  Figure  3.7  and  click  on  the  clip 
feature  to  select  it. 

4.  Select  File  >  Export  >  Export  Vector  Format.  From  the  pop-up  box, 
select  Shapefile. 

5.  A  dialog  box  will  appear  as  shown  in  Figure  3.16.  Click  Export  Areas.  A 
dialog  box  will  pop  up  in  which  to  type  in  a  file  name  to  save  the  shapefile. 
Make  sure  the  Generate  Projection  box  is  checked. 

6.  Click  on  the  Export  Bounds  tab  shown  in  Figure  3.17  and  click  on  the 

Crop  to  Selected  Area  Feature(s)  option.  Click  OK 
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Figure  3.15.  Changing  dispiay 
projection  to  UTM. 


Figure  3.16.  Diaiog  to  seiect  features 
forshapefiie. 


Shapefile  Export  Options 


File  Selection  j  Gridding  |  Export  Bounds  | 


Export  Areas 

Fie:  |M:\Ctel  XPS  Recover\Docum 
I”  Export  Lines 

Fie;  I  Select...  | 

I”  Export  Points 

Fie:  I  Select... 


The  shapefle  format  does  not  allow  the  mbdng  of  different 
entitiy  types  (.e.  points,  lines,  and  areas)  within  a  single 
shaF>e  file.  Thus,  you  need  to  specify  which  entity  types 
you  wish  to  export  and  what  file  to  export  them  to. 


Split  Expxxt  Based  on:  |  Do  Not  Split  Export  ^ 

Generate  projection  (P  FU)  file  for  each  entity  type .  Note 
R  that  some  software  ^.e.  Arc  Explorer)  cannot  handle 

shapefles  wth  a  PRJ  fie.  so  disable  if  you  have  trouble. 
r~  Generate  30  Features  Using  Loaded  Bevatkxi  Data 
W  Add  Feature  Type  (LAYER)  Attribute  to  DBF 
I”  Add  Feature  Map  Name  (MAP_NAME)  Attribute  to  DBF 
r  Add  Style  Attributes  to  DBF 
I”  Generate  Multi-Patch  Features  for  /Veas 


OK  I  Carxiel  |  Apply  |  Help 
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Figure  3.17.  Dialog  to  select  export 
bounds  for  shapefile. 


3.3.3  Delete  and  add  custom  buildings 

Some  of  the  building  details  have  changed  since  the  data  set  produced  for 
the  buildings  was  constructed  from  2009  imagery.  Some  buildings  that 
were  next  to  U.S.  Highway  75  were  deleted,  and  the  George  W.  Bush 
Presidential  Library  was  added,  as  shown  in  Figure  3.18. 

To  delete  features  in  Global  Mapper: 

1.  Select  the  building  using  the  Feature  Info  tool  (Figure  3.7). 

2 .  Click  the  Delete  button. 
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Figure  3.18.  Additional  edits  to  the  building  data  set. 


The  George  W.  Bush  Presidential  Library  was  added  to  the  building  data 
set  using  imagery  from  the  USDA  NAIP  from  2012,  as  well  as  images  from 
Google  (Figure  3.19)  for  reference.  To  add  the  custom  building,  perform 
the  following  steps  in  Global  Mapper: 

1.  Turn  the  imagery  layer  on  using  the  Overlay  Control  Center  as  shown 
in  Figure  3.5. 

2.  Select  the  Digitizer  Tool  and  click  the  Create  New  Area  Feature  as 
shown  in  Figure  3.20. 

3.  Click  on  points  to  trace  the  outline  of  the  building  as  shown  in  Figure  3.21. 
Right-click  on  the  last  point  to  end  the  feature  creation. 

4.  The  Modify  Feature  Info  dialog  box  will  display  as  shown  in  Figure  3.22. 
Click  on  Add. 

5.  The  Edit  Attribute/V alue  dialog  will  display  as  shown  in  Figure  3.23. 
Enter  the  attribute  Height  and  the  value  59.  The  value  of  59  ft  was 
considered  an  average  height  for  the  entire  building  from  information 
about  heights  contained  on  the  internet  and  inferred  from  Google  Earth. 
Click  on  OK. 
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Figure  3.19.  Image  of  Bush  Library  on  Google  Earth. 


Figure  3.20.  Create  new  area  feature. 

^  Ghsbat  Mapper  vl4.0  (b0^l2l2}  (64-biit]  -  REGISTERED  (Univers%Par^^l«jpg||^^^!lumerrto^prni^ 


File  Edit  View  T^Jol^  Terrain  Analysis  Search  GPS  Help 

e  n  tt 

%  ^  ^  -ft 

*ci  S  0‘  — 1  nl 

■Ki|4l4l 

<5.  K 1  .<■ 

Create  new  area  feature 


Select  digitizer  tool 


Figure  3.21.  Creating  new  area  feature  by 
tracing  photo. 
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Figure  3.22.  Name  the  added  feature. 

Modify  Feature  Info  ^ 


I  Bush  Library 


OK  I  Cancel  | 


Figure  3.23.  Add  an  attribute  to  feature. 


3.3.4  Trim  data  set  based  on  attributes 

It  was  necessary  to  trim  the  larger  building  data  set  into  smaller  data  sets 
that  satisfied  certain  criteria.  This  was  done  to  create  data  sets  that  only 
contained  buildings  of  a  certain  size  (i  m^)  and  height  (2  m).  This  was 
done  to  remove  smaller  features  such  as  dormers  and  ensure  that 
buildings  were  greater  than  a  specified  height.  Since  the  grid  spacing  of  the 
finite  difference  grid  was  approximately  1  m,  the  buildings  needed  to  be  at 
least  2  m  in  height  to  ensure  that  at  least  2  nodes  were  used  to  define  a 
building. 

In  ArcMap,  a  new  smaller  data  set  satisfying  certain  user-specified 
conditions  can  be  created  from  a  larger  data  set  using  the  following  steps: 

1.  Select  Selection  >  Select  by  Attributes  from  the  main  menu.  The 
dialog  box  shown  in  Figure  3.24  will  appear. 

2.  Select  the  Layer  that  is  to  be  used  to  perform  the  selection. 

3.  Select  the  Create  a  new  selection  option  for  the  Method. 
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4.  Create  the  selection  query  by  typing  the  requirements  into  the  text  box  as 
shown  in  Figure  3.24.  In  this  case,  all  buildings  are  selected  that  are 
greater  than  im^  in  area  and  greater  than  2  meters  (6.562  ft)  in  height.  The 
height  attribute  is  defined  in  feet. 

5.  Press  OK . 

6.  Right-click  on  the  layer  used  for  the  query  as  shown  in  Figure  3.25. 

7.  Select  Selection  >  Create  Layer  from  Selected  Features. 

8.  The  created  layer  will  appear  in  the  layer's  list  and  can  be  renamed  if 
desired. 

Once  a  new  layer  is  created,  it  may  also  be  used  to  further  refine  the  data. 
Therefore,  the  layer  created  from  the  above  procedure  may  be  further 
refined  to  segment  the  buildings  into  sets  that  are  less  than  or  equal  to  51 
ft  in  height  and  greater  than  51  ft  in  height.  This  height  was  selected  as  the 
height  that  best  segmented  the  buildings  into  residential  and  commercial 
buildings. 


Figure  3.24.  Select  by  attributes  to 
create  a  new  layer. 
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Figure  3.25.  Create  a  new  layer  from  the  selection. 
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3.3.5  Select  features  by  location 

Using  ArcMap,  data  sets  were  created  for  buildings  that  were  close  to  or 

away  from  U.S.  Highway  75.  To  create  the  data  set  that  contains  buildings 

within  200  m  of  U.S.  Highway  75,  perform  the  following  steps: 

1.  Select  Selection  >  Select  by  Location  from  the  main  menu.  The  dialog 
box  shown  in  Figure  3.26  will  appear. 

2.  Select  Select  Feature  From. 

3.  Select  the  Target  Layer  that  is  to  be  used  to  perform  the  selection.  This  is 
the  layer  from  which  the  objects  are  selected. 

4.  Select  the  Source  Layer  that  is  to  be  used  to  perform  the  selection.  This 
is  the  layer  that  objects  are  close  to. 

5.  Select  Target  Layer(s)  Features  are  Within  a  Distance  of  the 
Source  Layer  Feature. 

6.  Enter  200  m.  Press  OK  or  Apply. 

7.  As  shown  in  Figure  3.25,  right-click  on  the  layer  used  for  the  query  and 
select  Selection  >  Create  Layer  from  Selected  Features. 

8.  The  created  layer  will  appear  in  the  layer's  list  and  can  be  renamed  if 
desired. 
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Figure  3.26.  Select  buildings  (cyan  color)  within  200  m  of  U.S. 
Highway  75. 


To  create  the  data  set  that  contains  buildings  greater  than  200  m  from 
U.S.  Highway  75,  perform  the  following  steps: 

1.  As  shown  in  Figure  3.27,  right-click  on  the  layer  used  for  the  query  and 
select  Selection  >  Switch  Selection. 

2.  As  shown  in  Figure  3.25,  right-click  on  the  layer  used  for  the  query  and 
select  Selection  >  Create  Layer  from  Selected  Features. 

3.  The  created  layer  will  appear  in  the  layer's  list  and  can  be  renamed  if 
desired. 

3.3.6  Create  aggregated  building  sets 

The  effect  of  the  size  of  the  buildings  was  investigated  by  combing 
buildings  within  a  certain  distance  (aggregation  distance)  into  a  larger 
building  polygon.  The  aggregation  function  in  ArcMap  combined  smaller 
buildings  into  larger  ones  while  also  keeping  the  orthogonal  nature  of  the 
building  footprints  intact.  The  tools  in  ArcMap  made  these  computations 
easier  than  Global  Mapper.  The  aggregation  functions  could  be  performed 
in  Global  Mapper  by  creating  buffers,  but  the  ArcMap  functions  were 
much  quicker  and  gave  better  results. 
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Figure  3.27.  Switch  seiection  of  buiidings  (cyan  coior)  to  produce  buiidings 
away  from  U.S.  Highway  75. 


To  run  the  aggregation  function  in  ArcMap,  do  the  following: 

1.  Select  the  (ieoprocessing  >  ArcToolbox  option  from  the  main  menu. 

2.  Within  the  ArcToolbox,  select  Cartography  Tools  >  Generalization  > 
Aggregate  Polygons  as  shown  in  Figure  3.28. 

3.  From  the  dialog  shown  in  Figure  3.29,  select  the  Input  Features  to  use 
from  the  drop-down  menu. 

4.  Type  in  an  Output  Feature  Class  name. 

5.  Type  in  an  Aggregation  Distance  of  25  m. 

6.  Enter  a  Minimum  Area  of  o.  This  will  cause  the  algorithm  to  use  all 
buildings. 

7.  Enter  a  Minimum  Hole  Size  of  2500.  This  will  cause  the  algorithm  to 
discard  all  holes  below  this  minimum  size.  Holes  can  occur  in  the  larger 
aggregated  polygons  based  on  building  proximity. 

8.  Check  the  box  to  Preserve  Orthogonal  Shape.  This  will  cause  the 
algorithm  to  keep  the  outside  shape  of  the  building  footprints  when 
combining. 

9.  Press  OK  The  buildings  will  be  aggregated  as  shown  in  Figures  3.30  and 
3-31- 

Figure  3.30  shows  the  individual  buildings  overlaid  upon  the  aggregated 
polygons,  while  Figure  3.31  shows  just  the  aggregated  polygons.  This 
building  data  set  was  for  the  buildings  located  more  than  200  m  from 
U.S.  Highway  75.  The  same  aggregation  procedure  can  be  used  with  the 
building  data  set  containing  buildings  within  200  m  of  U.S.  Highway  75. 
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Figure  3.28.  Aggregate  Polygons  function  inside  ArcToolbox. 
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Figure  3.29.  Parameters  for  aggregating  buildings. 
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Figure  3.30.  Individual  (green)  and  aggregated  (purple)  buildings  located  away  from  U.S. 

Highway  75. 


Figure  3.31.  Aggregated  buildings  located  away  from  U.S.  Highway  75. 
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3.3.7  Compute  weighted  averages  of  heights  of  aggregated  buiidings 

Once  aggregated  building  data  sets  were  constructed,  an  appropriate 
height  attribute  still  needed  to  be  assigned  to  the  aggregated  polygons.  A 
weighted  average  of  the  height  was  computed  for  each  aggregated  polygon 
using  Equation  3.1: 


HtWtAvg 


E(AxHf) 

^4 


(3.1) 


where 

HtWtAvg=  weighted  average  of  heights  for  buildings  within  aggregated 
polygon 

A  =  area  of  individual  building  within  aggregated  polygon 
Ht  =  height  of  individual  building  within  aggregated  polygon 

To  compute  the  weighted  average  of  the  height  of  each  aggregated 

polygon: 

1.  Select  the  Geoprocessing  >  ArcToolbox  option  from  the  main  menu. 

2.  Within  the  ArcToolbox,  select  Analysis  Tools  >  Overlay  >  Spatial 
Join  as  shown  in  Figure  3.32. 

a.  From  the  dialog  shown  in  Figure  3.33,  select  the  Target 
Features  to  use  from  the  drop-down  menu. 

b.  Select  the  Join  Features  from  the  drop-down  menu. 

c.  Enter  the  Output  Feature  Class  to  save  the  query  to. 

d.  Select  JOIN_ONE_TO_MANY  for  the  Join  Operation. 

e.  Select  CONTAIN S  for  the  Match  Option. 

f.  Press  OK 

3.  Right-click  the  new  layer  created  and  select  Open  Attribute  Table.  The 
table  is  shown  in  Figure  3.34.  The  data  set  now  contains  each  individual 
building  with  its  corresponding  area  and  height  grouped  by  the 
TARGET_FID  which  defines  the  aggregated  polygon  for  the  individual 
buildings. 

4.  Click  the  left  icon  (Table  Options)  on  the  menu  of  the  Attribute  Table 
and  select  Add  Field.  This  will  be  used  to  compute  the  area  times  the 
height  of  each  building.  The  dialog  in  Figure  3.35  will  display. 

a.  Enter  AxHt  for  the  name  of  the  field. 
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b.  Select  Float  for  the  Type  of  field. 

c.  Set  the  field  properties  as  shown  in  Figure  3.35. 

d.  Press  OK  The  field  is  added  to  the  table  attributes  and  is 
initially  assigned  a  <Null>  value. 

5.  Right-click  the  newly  created  field  and  select  Field  Calculator  as  shown 
in  Figure  3.36.  The  Field  Calculator  will  display  as  shown  in  Figure  3.37. 

a.  Select  VB  Script. 

b.  Click  area  in  the  Fields  box. 

c.  Click  the  multiplication  button. 

d.  Click  Height  in  the  Fields  box. 

e.  Click  OK  This  will  compute  the  formula  entered  and  fill  in 
the  values  for  the  field  AxHt  as  shown  in  Figure  3.38. 

6.  A  summary  table  must  be  computed  to  sum  the  areas  and  the  areas  times 
the  heights  for  every  building  in  an  aggregated  polygon.  In  the  Attribute 
Table,  right-click  the  TARGET_FID  field  and  select  Summarize.  The 
Summarize  dialog  will  display  as  shown  in  Figure  3.39. 

a.  Select  TARGET_FID  from  the  Select  a  field  to 
summarize. 

b.  Check  the  Sum  box  for  both  the  area  and  the  height  as 
shown  in  Figure  3.39. 

c.  Under  Specify  output  table,  enter  a  name  for  the 
summary. 

d.  Click  OK. 

7.  The  summary  table  computed  must  be  joined  with  the  table  containing  the 
aggregated  buildings.  Right-click  the  layer  containing  the  aggregated 
buildings  in  the  Table  of  Contents. 

8.  Select  Joins  and  Relates  >  Joins.  The  Join  Data  dialog  shown  in 
Figure  3.40  will  display. 

a.  Select  Join  attributes  from  a  table. 

b.  Select  OB  JECTID  for  the  field  the  join  will  be  based  on. 

c.  Choose  the  summaiy  table  previously  computed  in  Step  6  as 
the  join  layer. 

d.  Select  OB  JECTID  as  the  field  to  base  the  join  on. 

e.  Select  Keep  all  records. 

f.  Press  OK.  The  attribute  table  shown  in  Figure  3.41  will  be 
created.  In  this  table,  the  sum  of  the  areas  and  the  sum  of 
the  areas  times  the  heights  has  been  created  for  each 
aggregated  polygon.  For  each  TARGET_FID  in  the  table 
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there  is  a  FREQUENCY  attribute  that  tells  how  many 
buildings  are  in  the  aggregated  polygon. 

9.  In  the  attribute  table  for  the  aggregated  polygons,  a  field  must  be  added  to 
compute  the  weighted  average  of  the  heights. 

a.  Add  a  field  to  the  table  as  described  in  Step  4  and  name  it 
WtAvgHt  as  shown  in  Figure  3.42. 

b.  Fill  in  the  rest  of  the  values  shown  in  Figure  3.42  and  click 

OK 

10.  Use  the  field  calculator  on  the  new  field  as  described  in  Step  5.  Fill  in  the 
values  shown  in  Figure  3.43. 

11.  The  results  of  the  field  calculator  showing  the  weighted  averages  of  the 
heights  for  the  individual  buildings  contained  in  an  aggregated  polygon  are 
shown  in  Figure  3.44. 

12.  Save  the  building  data  to  a  shapefile  as  described  in  Section  3.3.9. 


Figure  3.32.  Spatial  Join  function  in  ArcToolbox. 
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Figure  3.33.  Use  a  spatial  join  to  make  a  data  set  containing  building  height  information 

grouped  by  aggregated  polygon. 


Figure  3.34.  Result  of  spatial  Join  showing  heights  and  areas  of  buildings  grouped  by 

aggregated  polygon. 
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Figure  3.35.  Add  field  dialog. 


Figure  3.36.  Selection  of  the  field  calculator. 
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Figure  3.37.  Field  calculator  to  compute  area  times 
height  for  each  building. 
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Figure  3.38.  Result  of  calculations  using  field  calculator. 
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Figure  3.39.  Summarize  data  for  a  field  to  create  a 
new  table. 


Summarize  creates  a  new  table  containing  one  necond  for  each  unique  value 
of  the  selected  field,  along  with  statistics  summarizing  any  of  the  other  fields. 


1 .  Select  a  field  to  summarize: 

TARGET_FID 

2.  Choose  one  or  more  summary  statistics  to  be  included  in  the 
output  table; 


0  GM_TYPE 
0  ELEVATION 
0  Height 

EH  Minimum 
EH  Maximum 
EH  Average 
0  Sum 

EH  Standard  Deviation 
EH  Variance 
0  area 


S.  Specily  output  table: 


as  County\UniverBityPark\25C0m  area\Sum_0utpLrt_7.dbf 


O^iui^arize  on  the  sd^^ 


About  Summarizing  Data 


OK 


Cancel 


Figure  3.40.  Join  individuai  buiiding  data  with 
the  aggregated  buiiding  data. 
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Figure  3.41.  Result  of  join. 


Figure  3.42.  Add  a  field  for  the  weighted 
average  of  the  heights. 


ERDC  TR-15-5 


45 


Figure  3.43.  Use  the  field  calculator  to  compute  a 
weighted  average. 


Figure  3.44.  Table  showing  the  weighted  average  for  the  heights  of  the  buildings 
contained  in  aggregated  polygons. 
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3.3.8  Save  building  data  to  shapefiles  in  Global  Mapper 

After  clipping  operations  as  described  in  Section  3.2.2  or  addition  or 
deletion  of  buildings  as  described  in  Section  3.3.3,  the  data  set  can  be 
output  to  the  ESRI  shapefile  format  to  provide  one  file  that  contains  all 
editing  operations  for  later  use  with  ArcMap. 

The  output  of  the  data  is  the  same  as  described  in  Section  3.2.2.  The 
Export  Bounds  tab  as  shown  in  Figure  3.45  may  be  used  to  select  the 
bounds  of  the  data  to  export.  The  bounds  of  the  data  to  export  may  be 
specified  by  selecting: 

1.  Crop  to  Selected  Area  Feature(s)  which  will  crop  (clip)  the  data  to  the 
bounds  of  the  feature. 

2.  All  Loaded  Data  which  will  export  all  layers  turned  on  in  the  Overlay 
Control  Center. 

3.  All  Data  Visible  on  Screen  which  will  export  only  the  data  that  are 
visible  on  the  screen. 

Figure  3.45.  Shapefile  export  options  for 
Global  Mapper. 
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3.3.9  Save  building  data  to  shapefiles  in  ArcMap 

ArcMap  was  used  to  perform  editing  operations  on  the  building  data  set  in 
order  to  produce  several  data  sets  based  on  height,  area,  and  location  as 
explained  in  Sections  3.3.4  -  3.3.6.  To  save  the  data  sets  for  further 
processing  in  Global  Mapper,  the  data  can  be  saved  to  the  shapefile 
format.  To  save  data  in  ArcMap  to  the  shapefile  format; 

1.  Right-click  the  layer  of  interest  as  shown  in  Figure  3.46  and  select  Export 
Data. 

2.  The  dialog  in  Figure  3.47  will  appear.  Select  All  Features  and  type  a  file 
name  to  save  the  data  to.  Click  OK 

Figure  3.46.  Selecting  export  data  option  to  save  data. 
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Figure  3.47.  Export  data  dialog  for  shapefiles  in  ArcMap. 


3.3.10  View  buildings  in  3-D 

In  Global  Mapper,  the  building  data  may  be  viewed  in  g-D.  It  is  very  useful 
to  see  the  height  of  buildings  in  relation  to  the  terrain. 

The  building  data  from  CyberCity  3D  provided  the  height  attribute  of  the 
buildings  in  feet.  In  Global  Mapper,  the  attribute  used  to  select  the  elevation 
of  the  buildings  and  the  units  must  be  specified.  To  set  the  options: 

1.  In  the  Overlay  Control  Center,  click  on  a  layer  containing  the  building 
data.  Turn  on  other  layers  of  interest  to  display  in  the  3-D  view  such  as 
highways  or  orthophotography.  Turn  on  the  topography  layer  to  display 
the  buildings  on  the  actual  terrain;  otherwise,  the  buildings  will  be 
displayed  on  a  flat  surface. 

2.  Click  the  Options  button.  The  Vector  Options  dialog  will  display  as 
shown  in  Figure  3.48. 

a.  Set  the  Get  Elevation  from  Attribute  Value  to 
HEIGHT. 

b.  Set  the  Elevation  Units  for  Unspecified  Elevation 
Values  to  FEET. 

3.  Select  the  3D  View  icon  on  the  main  menu  to  display  the  3-D  view  as 
shown  in  Figure  3.49.  The  3D  view  will  display  as  shown  in  Figure  3.50. 

4.  Select  the  Change  display  properties  from  the  main  menu  of  the  3-D 
view  as  shown  in  Figure  3.50.  The  dialog  to  change  display  properties  will 
appear  as  shown  in  Figure  3.51. 
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a.  Under  3D  Vector  Display  Options,  check  all  boxes. 

b.  Since  the  height  of  the  buildings  is  specified  as  a  height  and 
not  an  elevation,  the  box  Treat  3D  Vector  Elevations  as 
Relative  to  the  Terrain  Surface  must  be  checked. 

c.  Checking  the  Extrude  3D  Areas  to  Surface  causes  the 
walls  of  the  buildings  to  be  extruded  to  the  surface  of  the 
terrain. 


Figure  3.48.  Vector  options  diaiog  to  set 
eievation  data  for  buiidings. 


Figure  3.49.  Seiect  3D  view. 


^  Global  Mapper  vl4.0  (b[>91212)  [64-bit]  -  REGISTERED  (UniversityPark-25CIOarea-for  report.gmw) 
File  Edit  View  Tools  Terrain  Analysis  Search  GPS  Help 
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Figure  3.50.  Change  display  properties  of  3D  view. 


Figure  3.51. 3D  view  properties. 


3.3.11  Create  files  for  Urban  Modeler 

A  piece  of  custom  software,  termed  the  Urban  Modeler,  was  used  to 
provide  further  processing  of  data  obtained  from  Global  Mapper.  The 
vector  information  concerning  buildings  can  be  output  from  Global 
Mapper  in  several  vector  file  formats.  The  output  of  the  data  in  the 
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shapefile  format  is  described  in  Section  3.3.8.  The  ASCII  vector  file  format 
was  chosen  for  ease  of  use  with  the  Urban  Modeler  software. 

To  output  data  in  the  ASCII  vector  file  format  from  Global  Mapper: 

1.  The  data  should  be  exported  in  the  UTM  coordinate  system  for  use  with 
Urban  Modeler.  If  the  coordinate  projection  needs  to  be  changed,  follow 
the  instructions  in  Section  3.2.3  and  select  the  UTM  coordinate  projection. 

2.  Select  the  Overlay  Control  Center  button  from  the  main  toolbar. 

3.  Turn  on  the  building  data  layer(s)  by  checking  the  appropriate  layers. 

4.  From  the  main  menu,  select  File  >  Export  >  Export  Vector  Format. 

5.  From  the  dialog  shown  in  Figure  3.52,  select  Simple  ASCII  Text  File. 

6.  The  ASCII  Export  Options  dialog  will  appear  as  shown  in  Figure  3.53. 

a.  Under  Coordinate  Separator,  click  Comma. 

b.  Under  Feature  Separator,  click  Custom  and  then  enter 
\nBuilding.  This  causes  a  new  line  and  the  word  Building 
to  be  written  before  each  data  record.  Urban  Modeler  uses 
sets  of  key  words  to  recognize  types  of  data. 

c.  Check  Export  Elevation  for  Each  Vertex. 

d.  Check  Include  Feature  Attributes  Before  Coordinate 
Data. 

7.  The  default  options  on  the  other  tabs  may  be  left  as  the  default. 

8 .  Click  OK  and  input  a  file  name. 

Figure  3.52.  Export  file  options. 
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Figure  3.53.  ASCII  export  options  for 
buildings. 


The  building  data  will  be  output  in  one  file.  Urban  Modeler  will  read  the 
data  and  perform  additional  processing  as  discussed  in  Chapter  4. 

The  units  of  the  UTM  coordinate  system  are  meters,  whereas,  the  units  of 
the  height  attribute  for  the  buildings  defined  by  CyberCity  3D  are  in  feet. 
Global  Mapper  does  not  alter  attribute  units;  therefore,  the  height 
attribute  will  remain  in  feet  even  though  the  building  footprint  data  are 
defined  in  meters. 

In  Figure  3.53,  the  option  to  export  the  elevation  for  each  vertex  may  be 
selected.  This  will  cause  the  height  attribute  to  be  written  as  the  z- 
coordinate  and  it  will  be  converted  to  meters  automatically. 

Global  Mapper  also  allows  the  z-coordinate  to  be  set  from  the  topography. 
This  option  is  not  used  for  buildings,  but  is  for  bridges  as  discussed  in 
Section  3.5.  If  the  z-coordinate  of  the  vertices  of  an  object  is  set  using  the 
topographic  elevation,  the  units  of  the  coordinates  will  be  converted  to 
meters. 
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3.4  Model  paved  areas 

3.4.1  Read  paved  areas  into  Global  Mapper 

Global  Mapper  was  used  to  read  2-D  planimetric  data  from  NCTCOG 
defining  the  paved  areas  of  the  urban  area  shown  in  Figure  3.54  in  the 
ESRI  shapefile  format.  The  coordinates  of  the  lower-left  corner  are  shown 
in  UTM  coordinates  in  Figure  3.54. 


Figure  3.54.  Paved  areas  of  analysis  region. 


Data  are  read  into  Global  Mapper  by  opening  the  specific  file  using  the 
File  menu  option  or  dragging  the  file  from  Windows  Explorer.  If  the  file 
contains  a  definition  of  the  projected  coordinate  system  used,  Global 
Mapper  will  just  automatically  read  and  display  the  data.  If  a  projected 
coordinate  system  is  not  found,  a  dialog  box  will  display  asking  the  user  to 
verify  the  appropriate  projected  coordinate  system  as  shown  in  Figure  3.1. 
The  user  would  select  the  appropriate  projected  coordinate  system  that  the 
data  set  uses  and  select  the  OK  button. 

3.4.2  Clip  paved  areas  to  desired  area 

The  initial  planimetric  data  set  obtained  from  NCTCOG  was  much  larger 
than  the  analysis  region.  Therefore,  the  data  set  was  clipped  at  the 
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boundaries  of  the  analysis  region  as  shown  by  the  green  outline  in 
Figure  3.54.  The  clipping  procedure  is  the  same  as  described  for  buildings  in 
Section  3.3.2.  The  result  of  the  clipping  operation  was  an  ESRI  shapefile 
that  contained  only  the  paved  areas  within  the  green  outline  shown  in 
Figure  3.54. 

3.4.3  Break  larger  paved  areas  into  smaller  areas 

When  processing  area  features  such  as  the  paved  areas,  it  was  found  that 
the  MATLAB  preprocessing  algorithms  performed  better  if  the  data  were 
split  into  smaller  areas.  This  seemed  to  help  with  the  speed  of  processing 
plus  alleviated  any  memory  problems  arising  from  trying  to  process  large 
areas. 

Breaking  larger  area  features  into  smaller  area  features  is  accomplished  in 
Global  Mapper  when  saving  the  vector  data  to  a  vector  output  file  (i.e.,  a 
shapefile  or  ASCII  vector  file)  as  described  in  Section  3.4.5. 

3.4.4  Detect  islands  in  data 

Islands  represent  open  spaces  in  a  larger  area.  That  is,  a  building  may  exist 
with  an  interior  courtyard,  which  can  be  modeled  using  an  island.  Paved 
areas,  such  as  highways,  may  have  interior  open  spaces  that  represent  the 
medians.  The  network  of  streets  may  be  represented  as  large  solid  areas 
with  islands  representing  the  city  blocks.  When  data  containing  islands  are 
displayed  in  GIS  programs,  the  islands  will  cause  the  large  solid  areas  to 
appear  to  have  open  spaces  in  them. 

For  the  MATLAB  preprocessing  routines  to  properly  handle  the  islands, 
the  islands  must  be  detected  and  written  out  separately.  This  is 
accomplished  with  the  Urban  Modeler  program.  No  special  processing  is 
needed  within  Global  Mapper  to  account  for  the  islands.  To  detect  the 
presence  of  islands  in  Global  Mapper: 

1.  Select  the  Overlay  Control  Center  button  from  the  main  toolbar. 

2.  Click  the  layer  containing  the  paved-area  data  to  view  it  on  the  screen. 

3.  Right-click  the  layer  containing  the  paved-area  data  in  the  Overlay 
Control  Center. 

4.  From  the  pop-up  menu,  choose  Select  All  Features  in  Selected 
Layer(s)  with  Digitizer  Tool.  The  points  defining  the  data  will 
highlight. 
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5.  Right-click  on  the  highlighted  data  on  screen. 

6.  Select  Advanced  Selection  Options  >  Select  Island  Areas  in 
Selected  Area  Feature. 

7.  The  display  will  highlight  all  islands  in  red  as  shown  in  Figure  3.55. 


Figure  3.55.  Island  areas  contained  in  the  paved  area  data  highlighted  in  red. 


3.4.5  Create  files  for  Urban  Modeler 

A  piece  of  custom  software,  termed  the  Urban  Modeler,  was  used  to 
provide  further  processing  of  data  obtained  from  Global  Mapper.  The 
vector  information  concerning  paved  areas  can  be  output  from  Global 
Mapper  in  several  vector  file  formats.  The  output  of  the  data  in  the 
shapefile  format  is  described  in  Section  3.3.8.  The  ASCII  vector  file  format 
was  chosen  for  ease  of  use  with  the  Urban  Modeler  software. 

To  output  data  in  the  ASCII  vector  file  format  from  Global  Mapper: 

1.  The  data  should  be  exported  in  the  UTM  coordinate  system  for  use  with 
Urban  Modeler  and  the  MATLAB  preprocessing  routines.  If  the  coordinate 
projection  needs  to  be  changed,  follow  the  instructions  in  Section  3.2.3  and 
select  the  UTM  coordinate  projection. 

2.  Select  the  Overlay  Control  Center  button  from  the  main  toolbar. 

3.  Turn  on  the  paved  area  data  layer(s)  by  checking  the  appropriate  layers. 
Also  turn  on  the  clipping  area. 
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4.  Select  the  Feature  Info  Tool  as  shown  in  Figure  3.7  and  click  on  the  clip 
area  to  select  it. 

5.  From  the  main  menu,  select  File  >  Export  >  Export  Vector  Format. 

6.  From  the  dialog  shown  in  Figure  3.52,  select  Simple  ASCII  Text  File. 

7.  The  ASCII  Export  Options  dialog  will  appear  as  shown  in  Figure  3.56. 

a.  Under  Coordinate  Separator,  click  Comma. 

b.  Under  Feature  Separator,  click  Custom  and  then  enter 
\nPaved  Area.  This  causes  a  newline  and  the  words 
Paved  Area  to  be  written  before  each  data  record.  Urban 
Modeler  uses  sets  of  key  words  to  recognize  types  of  data. 

c.  Check  Export  Elevation  for  Each  Vertex. 

d.  Check  Include  Feature  Attributes  Before  Coordinate 
Data. 

8.  Click  on  the  Gridding  tab  shown  in  Figure  3.56.  The  tab  shown  in  Figure 
3.57  will  display. 

9.  Click  Specified  Number  of  Rows  and  Columns.  Input  the  number  of 
rows  and  columns  desired.  As  shown  in  Figure  3.57,  enter  5  rows  and  5 
columns.  The  data  will  be  segmented  into  a  grid  with  the  data  split  along 
the  grid  boundaries.  For  a  5  by  5  grid,  25  files  will  be  written. 

10.  The  grid  naming  convention  can  be  left  as  shown  in  Figure  3.57. 

11.  Click  on  the  Export  Bounds  tab  shown  in  Figure  3.58  and  click  on  the 
Crop  to  Selected  Area  Feature(s)  option. 

12.  Click  OK  and  input  a  file  name. 

Figure  3.56.  ASCII  export  options  for 
paved  areas. 
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Figure  3.57.  Gridding  option  for 
paved  areas. 


Figure  3.58.  Export  bounds  options 
for  paved  areas. 
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3.5  Model  bridges 

3.5.1  Select  bridges  from  paved  areas 

Bridges  were  contained  within  the  planimetric  data  obtained  from 

NCTCOG.  They  were  denoted  with  an  attribute  called  FEATURE,  which 

contained  the  values  of  BRIDGE  or  PEDESTRIAN  BRIDGE.  To  extract  the 

bridge  information  using  Global  Mapper: 

1.  Select  the  Overlay  Control  Center  button  from  the  main  toolbar. 

2.  Turn  on  the  paved  area  data  layer(s)  by  checking  the  appropriate  layers. 

3.  Click  the  Select  by  Attribute/Name/Description  button  on  the  main 
toolbar,  as  shown  in  Figure  3.59. 

4.  Sort  the  features  by  clicking  on  the  FEATURE  tab  as  shown  in 
Figure  3.60.  This  will  group  all  the  bridges  together. 

5.  Select  all  bridges  by  clicking  on  the  first  BRIDGE  object.  Scroll  down  to 
the  last  BRIDGE  object  and  while  holding  down  the  SHIFT  key,  click  the 
BRIDGE  object.  You  may  also  scroll  down  to  the  PEDESTRIAN 
BRIDGE  objects  and  individually  select  a  bridge  by  holding  down  the 
CTRL  key  and  clicking  on  each  object. 

6.  Right-click  the  selected  BRIDGE  objects  and  select  Copy  the  Selected 
Objects  to  the  Clipboard. 

7.  Click  on  Edit  >  Paste  Features  From  Clipboard.  The  dialog  shown  in 
Figure  3.61  will  display. 

8.  Select  Paste  to  New  Layer.  This  will  create  a  layer  containing  just  the 
bridges. 

9.  Right-click  the  new  bridge  layer  in  the  Overlay  Control  Center  and 
select  Select  All  Features  in  Selection  Layer(s)  with  Digitizer 
Tool. 

10.  Right-click  on  a  selected  bridge  and  select  Attribute  Functions  > 
Apply  Elevations  from  Terrain  Layers  to  Selected  Feature(s). 
This  adds  a  z-coordinate  to  the  vertices  of  the  bridge  object  based  on  the 
elevation  of  the  underlying  terrain.  These  elevations  will  be  used  in  Urban 
Modeler  to  compute  the  placement  of  the  bridge. 

Figure  3.59.  Select  features  by  attribute. 

IB  Global  Msporr  vlU  ;M-bit]  -  REGISHHED  for 


Select  by  Attribute/Name/Description 
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Figure  3.60.  Selecting  bridges  based  on  attribute  value. 


Figure  3.61.  Copy  objects  to  new  layer. 


3.5.2  Flatten  terrain  under  bridges 

The  LIDAR  data  from  TNRIS  resulted  in  raised  areas  of  terrain  under  the 
bridges,  which  should  be  flattened  (or  lowered)  to  the  elevation  of  the 
roadway  as  shown  in  Figure  3.62.  To  flatten  the  terrain  under  a  bridge: 


1.  Turn  the  imagery  layer  on  using  the  Overlay  Control  Center  as  shown 
in  Figure  3.5. 

2.  Select  the  Digitizer  Tool  and  click  the  Create  New  Area  Feature  as 
shown  in  Figure  3.20. 
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3.  For  each  bridge,  click  on  points  to  trace  the  outline  of  the  area  under  the 
bridge  to  be  flattened  as  shown  in  Figure  3.63.  Right-click  on  the  last  point 
to  end  the  feature  creation. 

4.  The  Modify  Feature  Info  dialog  box  will  display  as  shown  in 
Figure  3.64.  Click  OK 

5.  Continue  to  create  areas  defining  the  regions  to  flatten  beneath  each 
bridge  feature. 

6.  When  all  areas  are  defined,  click  and  drag  the  mouse  to  create  a  selection 
box  around  the  new  areas. 

7.  Right-click  on  one  of  the  selected  areas  and  select  Attribute  Functions 
>  Apply  Elevations  from  Terrain  Layers  to  Selected  Feature(s). 
This  adds  a  z-coordinate  to  the  vertices  of  the  areas  used  to  flatten  the 
terrain. 

8.  Double-click  on  an  individual  area  to  select  it.  The  dialog  shown  in 
Figure  3.65  will  appear. 

9.  Click  on  the  Vertices  button.  The  dialog  in  Figure  3.66  will  appear. 

10.  Click  on  the  first  vertex.  Move  down  to  the  last  vertex  and  press  the 
SHIFT  key  and  click  on  the  vertex.  This  will  select  all  vertices. 

11.  Click  the  Edit  Elevation  button.  The  dialog  in  Figure  3.67  will  appear. 

12.  Enter  an  elevation  and  click  OK.  All  elevation  values  will  change  to  the 
new  value. 

13.  Select  all  area  features  as  described  in  Step  6. 

14.  Right-click  an  area  feature  and  select  Advanced  Feature  Creation 
Options  >  TERRAIN  -  Create/Flatten  Terrain  from  Selected 
Area  Feature(s). 

15.  The  Elevation  Grid  Creation  Options  dialog  will  display  as  shown  in 
Figure  3.68. 

16.  Check  the  Flatten  3D  Area  Features.  Click  OK. 

17.  Figures  3.69  and  3.70  show  the  result  of  flattening  the  terrain  under  the 
bridge. 
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Figure  3.62.  Raised  terrain  under  bridge. 


Figure  3.63.  Create  area  feature  encompassing  area  to  flatten. 
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Figure  3.64.  Dialog  for  newly  created  area  feature. 


Figure  3.65.  Dialog  for  selected  area. 
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Figure  3.66.  Selection  of  vertices  to  edit. 


Figure  3.67.  Editing  elevation  of 
all  vertices. 


Figure  3.68.  Elevation  grid  options. 


Grid  Optioos  jGriddng  |  Grid  Bounds  | 


Vertical  Unts  |  METERS 
-Grid  Spacing  - 


(•  Automaticaly  Determine  Optimal  Grid  Spacing 
f'  Manualy  Specify  the  Grid  Spacing  to  Use 
X-aods:  [o  meters 

Y-axis;  [o  meters 

If  you  wish  to  change  the  ground  units  that  the  resolution  is 
specified  h.  you  need  to  change  the  current  projection  on  the 
Projection  tab  of  the  Configurabon  dialog. 


pBevabon  Grid  "No  Data"  Distance  Criteria - 

This  setting  controls  how  farfrom  a  known  data  point  that  an 
elevation  grid  cel  has  to  be  before  it  is  considered  invalid.  The 
default  setting  assumes  al  grid  points  are  valid.  Lower  values 
make  the  valid  grid  stay  tighter  around  known  data  points. 


Tight 


Loose 


R  Use  3D  Line  Features  as  Constraints  f  .e.  Breakiines) 

W  Ratten  3D  Area  Features 

r~  Taper  3D  Area  Features  Using  Curve  Value;  fl 

r~  Ignore  Zero  Bevations 

l~  Save  Triangulation  Network  (TIN)  as  a  Vector  Layer 
r~  Heights  Relative  to  Ground  (Using  Loaded  Grid  Layers) 

[~  Fi  Entire  Bounding  Box  Instead  of  Just  Inside  Convex  Hul 


1 

OK 

1  Cancel  | 

Apply  1 

Help 

ERDC  TR-15-5 


64 


Figure  3.69.  Flattened  terrain  under  bridge. 


Figure  3.70.  3-D  view  of  flattened  area. 
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3.5.3  Save  terrain  with  fiattened  areas  to  GeoTiFF 

The  terrain  with  the  flattened  portions  may  be  selected  and  output  as  a 
32-bit  GeoTIFF  file  as  described  in  Section  3.2.3. 

3.5.4  Create  fiies  for  Urban  Modeier 

A  piece  of  custom  software,  termed  the  Urban  Modeler,  was  used  to 
provide  further  processing  of  data  obtained  from  Global  Mapper.  The 
vector  information  concerning  bridge  features  can  be  output  from  Global 
Mapper  in  several  vector  file  formats.  The  output  of  the  data  in  the 
shapefile  format  is  described  in  Section  3.3.8.  The  ASCII  vector  file  format 
was  chosen  for  ease  of  use  with  the  Urban  Modeler  software. 

To  output  data  in  the  ASCII  vector  file  format  from  Global  Mapper: 

1.  The  data  should  be  exported  in  the  UTM  coordinate  system  for  use  with 
Urban  Modeler  and  the  MATLAB  preprocessing  routines.  If  the  coordinate 
projection  needs  to  be  changed,  follow  the  instructions  in  Section  3.2.3  and 
select  the  UTM  coordinate  projection. 

2 .  Select  the  Overlay  Control  Center  button  from  the  main  toolbar. 

3.  Turn  on  the  bridge  data  layer(s)  by  checking  the  appropriate  layers.  Also 
turn  on  the  clipping  area. 

4.  Select  the  Feature  Info  Tool  as  shown  in  Figure  3.7  and  click  on  the  clip 
area  to  select  it. 

5.  From  the  main  menu,  select  File  >  Export  >  Export  Vector  Format. 

6.  From  the  dialog  shown  in  Figure  3.52,  select  Simple  ASCII  Text  FUe. 

7.  The  ASCII  Export  Options  dialog  will  appear  as  shown  in  Figure  3.71. 

a.  Under  Coordinate  Separator,  click  Comma. 

b.  Under  Feature  Separator,  click  Custom  and  then  enter 
\nBridge.  This  causes  a  newline  and  the  words  Bridge  to 
be  written  before  each  data  record.  Urban  Modeler  uses  sets 
of  key  words  to  recognize  types  of  data. 

c.  Check  Export  Elevation  for  Each  Vertex. 

d.  Check  Include  Feature  Attributes  Before  Coordinate 
Data. 

e.  Check  Export  Elevation  for  Each  Vertex. 

f.  Check  Include  Feature  Attributes  Before  Coordinate 
Data. 

8.  The  default  options  on  the  other  tabs  may  be  left  as  the  default. 

9.  Click  OK  and  input  a  file  name. 
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Figure  3.71.  ASCII  export  options  for 
bridges. 


3.6  Model  other  area  features 

The  processing  discussed  in  this  chapter  was  specific  to  buildings,  paved 
areas,  and  bridges.  The  information  given  may  be  used  to  process  other 
types  of  area  features,  such  as,  water  features  or  forested  areas.  These 
features  were  not  specifically  modeled  in  the  study  described  in  this 
report. 
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4  Run  Urban  Modeler  to  Produce  CSV  Files 

A  piece  of  custom  software  called  the  Urban  Modeler  was  used  to  read  in 
the  sets  of  ASCII  vector  data  that  were  output  by  Global  Mapper  as 
described  in  Chapter  3.  The  processing  of  this  data  corresponds  to  Step  4 
in  the  process  map  of  Figure  1.4.  Urban  Modeler  was  used  to: 

•  Convert  the  units  of  vertical  measurement  of  attributes  to  a  common 
unit  of  meters. 

•  Check  if  islands  exist  in  the  data  and  write  the  island  information  out 
separately  from  the  feature  information. 

•  Break  large  sets  of  data  (e.g.,  thousands  of  buildings)  into  many 
smaller  sets  of  data  to  improve  processing  later  in  MATLAB. 

•  Process  bridge  information  to  compute  the  top  and  bottom  elevation  of 
the  bridge  features. 

•  Eliminate  objects  that  have  a  small  area  (less  than  1  m^)  or  height  (less 
than  2  m).  ArcMap  was  used  to  provide  this  same  functionality. 

•  Eliminate  islands  with  a  small  (less  than  1  m^)  or  zero  area. 

•  Produce  data  files  for  buildings  and  area  features  (i.e.,  paved  features 
and  bridges). 

•  Write  GIS  information  out  to  ASCII  comma  separated  variable  (CSV) 
files  for  use  with  Microsoft  Excel. 

Islands  represent  open  spaces  in  a  larger  area.  That  is,  a  building  may  exist 
with  an  interior  courtyard,  which  can  be  modeled  using  an  island.  Paved 
areas,  such  as  highways,  may  have  interior  open  spaces  that  represent  the 
medians.  The  network  of  streets  may  be  represented  as  large  solid  areas 
with  islands  representing  the  city  blocks.  The  islands  must  be  written  out 
separately  to  properly  account  for  the  material  properties  of  the  objects. 

The  CSV  files  produced  by  the  Urban  Modeler  program  define  the 
structures  (buildings  and  paved  areas)  that  sit  on  the  topography.  These 
files  are  used  in  Microsoft  Excel  to  define  the  object  geometry  and  material 
properties  for  use  with  the  MATLAB  processing  routines.  The  use  of  the 
Excel  spreadsheets  is  described  in  Chapter  5. 
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4.1  Produce  CSV  files  for  buildings 

The  building  information  produced  by  Global  Mapper  is  input  into  the 

Urban  Modeler  program  and  processed  as  follows: 

1.  Select  File  >  Import  >  Buildings  to  select  an  ASCII  data  file  from 
Global  Mapper  containing  building  information  as  shown  in  Figure  4.1. 

2.  The  dialog  in  Figure  4.2  will  display.  Click  OK. 

3.  Repeat  Steps  1  and  2  to  continue  to  read  in  building  files  until  all  building 
data  has  been  entered.  The  complete  building  data  set  is  shown  in 
Figure  4.3. 

4.  Select  FD  Grid  >  Save  Building  Footprints  to  CSV.  This  saves  the 
footprint  of  the  building  features  in  a  format  necessary  for  the  MATLAB 
preprocessing  routines. 

5.  Input  a  file  name. 

The  buildings  are  saved  to  a  CSV  file  for  later  use  in  Microsoft  Excel. 

Multiple  files  are  written  depending  upon  the  number  of  building  objects. 

Files  are  written  containing  3000  buildings  per  file.  The  file  name  contains 

the  name  of  the  CSV  file  entered  plus  an  additional  numeric  designator 

identifying  the  number  of  the  output  file  (e.g.,  buildings_away4.csv). 


Figure  4.1.  Reading  buiiding  information  into  the  Urban  Modeier  program. 


t. 
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Figure  4.2.  Defining  method  used  to 
denote  height  of  building  in  input  file. 


Figure  4.3.  Complete  building  data  set. 


4.2  Produce  CSV  files  for  areal  features 

The  area  feature  information  for  both  paved  area  and  bridge  features 
output  by  Global  Mapper  are  processed  by  Urban  Modeler  to  produce  files 
necessary  for  the  MATLAB  preprocessing  routines.  The  procedure  to 
produce  the  files  is  as  follows: 

1.  Select  File  >  Import  >  Area  Features  to  select  an  ASCII  data  file  from 
Global  Mapper  containing  the  paved  area  or  bridge  information.  Urban 
Modeler  will  ask  if  multiple  files  are  to  be  read  in  based  on  a  base  name. 
This  can  be  used  to  read  in  multiple  files  with  a  common  base  name.  For 
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example,  the  paved  area  information  consisted  of  25  data  sets  with  names 
ending  in  the  letters  A  through  E  and  numbers  1  through  5.  Therefore, 
there  were  five  files  with  the  letter  A,  five  files  with  the  letter  B,  and  so  on. 
An  example  file  name  would  be  PavedAreas_A4.txt.  To  read  in  all  files 
containing  the  letter  A,  the  user  would  enter  PavedAreas_A.txt  when 
prompted  for  a  file  name. 

2.  Repeat  Step  1  to  continue  to  read  in  paved  area  or  bridge  files  until  all  data 
have  been  entered.  The  user  can  read  in  a  single  input  file  and  process  the 
data  or  read  in  multiple  input  files  and  process  many  sets  of  data  at  once. 

3.  The  complete  paved  area  data  set  is  shown  in  Figure  4.4.  The  yellow 
shaded  areas  are  the  islands.  The  bridge  data  are  shown  in  Figure  4.5.  The 
paved  area  and  bridge  data  sets  should  be  processed  separately  so  that  the 
resulting  data  can  be  written  to  separate  files.  This  is  needed  because  the 
file  structure  is  not  exactly  the  same  for  both  sets  of  data. 

4.  Select  FD  Grid  >  Save  Area  Sets  to  CSV.  This  saves  the  footprint  of  the 
area  features  in  a  format  necessary  for  the  MATLAB  preprocessing 
routines.  The  program  will  detect  if  islands  are  present  and  write  out  two 
files.  One  file  will  contain  the  area  feature  definitions  while  the  other  file 
with  -islands  appended  to  the  file  name  will  contain  the  island 
definitions.  Urban  Modeler  will  ask  if  one  output  file  or  multiple  output 
files  are  to  be  written.  If  multiple  output  files  are  selected,  area  features 
will  be  written  to  separate  output  files  with  names  based  on  the  input  file 
names  read.  For  example,  if  five  input  files  had  been  read  in  with  the  name 
PavedAreas-A.txt  as  discussed  in  Step  1,  the  output  files  would  be 
named  PavedAreas_A.csv  and  PavedAreas_A-islands.csv. 

5.  Input  a  file  name  to  save  the  processed  data  sets. 

All  area  features  representing  roads,  streets,  highways,  parking  lots,  and 
medians  are  assigned  a  height  of  -2  m  when  processed.  The  negative 
height  causes  the  area  feature  to  be  placed  below  the  topography.  This  is 
explained  in  greater  detail  in  Chapter  5. 

Bridges  features  are  defined  by  a  perimeter.  The  points  defining  the 
perimeter  of  the  bridge  have  z-coordinates  taken  from  the  topography  and 
assigned  by  Global  Mapper  as  discussed  in  Section  3.5.1.  Urban  Modeler 
uses  these  z-coordinates  to  compute  a  minimum  elevation  for  the  bridge. 
The  thickness  of  the  bridge  is  set  to  2  m. 
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Figure  4.4.  Paved  area  data  set  in  Urban  Modeler. 
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Figure  4.5.  Bridge  information  displayed  in  Urban  Modeler. 
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5  Use  of  Microsoft  Excel  to  Describe  Data 

As  described  in  Chapter  4,  CSV  files  produced  by  the  Urban  Modeler 
program  are  used  to  aid  in  the  definition  of  the  urban  area.  The  CSV  files 
are  inserted  into  various  Excel  spreadsheets  to  aid  in  defining  auxiliary 
structures  such  as  buildings,  bridges,  and  paved  areas.  The  Excel 
spreadsheets  are  also  used  to  define  the  grid,  air  layers,  material 
properties,  and  topography.  The  spreadsheets  are  used  by  the  MATLAB 
preprocessing  routines  to  construct  the  data  files  necessary  to  run  the 
PSTOP3D  program.  This  chapter  explains  the  setup  of  the  spreadsheets 
and  corresponds  to  Step  5  of  the  process  flowchart  in  Figure  1.4. 

5.1  Open  CSV  files  in  Excel 

Information  contained  within  CSV  files  can  easily  be  put  into  an  Excel 
spreadsheet  by  opening  the  CSV  file  directly  into  Excel  and  copying  the 
information  into  the  appropriate  spreadsheet.  This  can  be  done  as  follows: 

1.  Open  the  master  Excel  workbook  containing  the  spreadsheet  templates. 

2.  Double-click  on  the  CSV  file.  The  file  will  open  in  Excel. 

3.  Click  on  the  first  line  of  data. 

4.  Press  CTRL-SHIFT-END  to  select  all  data. 

5.  Press  CTRL-C  to  copy  data. 

6.  Go  to  the  specific  spreadsheet  required  for  the  data  (e.g.,  a  building 
template)  and  press  CTRL-V  to  paste  the  data  into  the  sheet. 

7.  Repeat  the  above  procedure  for  each  CSV  file. 

5.2  General  description  of  use  of  Excel  spreadsheets 

Excel  spreadsheets  are  used  to  setup  the  definition  of  the  data  for  the 
acoustic  analysis.  The  spreadsheets  are  contained  in  an  Excel  workbook. 
There  are  several  types  of  spreadsheets  that  control  the  setup  of  the 
acoustic  analysis.  For  the  acoustic  analyses,  the  Excel  workbook  consists  of 
spreadsheets  to  define  the: 

1.  Problem  setup. 

2.  Size  of  the  grid  used  in  the  analysis. 

3.  Topography. 
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4.  Air  layers  and  type  of  material  properties  (e.g.,  poroacoustic)  used  in  the 
analysis. 

5.  Auxiliary  structures  (e.g.,  buildings,  paved  areas,  and  bridges). 

Some  of  the  spreadsheets  (e.g.,  the  spreadsheet  defining  buildings)  will 
consist  of  a  top  or  header  section  that  will  not  change  and  a  lower  section 
that  will  consist  of  the  information  contained  in  the  CSV  files  generated  by 
the  Urban  Modeler  program. 

The  spreadsheets  consist  of  various  data  items  designated  by  command 
words  with  an  associated  value  of  the  data  item.  The  command  words  for 
the  data  items  should  be  left  as  shown  in  the  spreadsheets. 

Each  of  these  spreadsheets  will  be  discussed  in  the  following  sections 
explaining  the  data  items  and  options  available  for  an  acoustic  analysis. 
The  data  items  and  options  that  are  explained  are  the  ones  used  in  this 
study  for  acoustic  analyses.  More  options  exist  than  will  be  covered  in  this 
discussion;  these  options,  for  example,  are  for  seismic  analysis.  Some 
options  may  be  left  blank  as  shown  in  the  following  sections. 

5.3  Main  spreadsheet  to  describe  probiem  setup 

The  spreadsheet  shown  in  Figure  5.1  contains  information  that  controls 
the  naming  of  data  files  and  the  names  of  the  spreadsheets  used  to  define 
the  problem. 

The  data  items  in  the  spreadsheet  are  defined  as  follows: 

1.  options.model.type  should  be  set  to  acoustic  to  compute  particle 
velocities  and  pressure  fields. 

2.  id.model  is  a  designator  for  the  model  being  run. 

3.  id.grid  is  the  name  of  the  grid. 

4.  id.topo  is  the  name  of  the  topography. 

5.  id.med  is  the  name  of  the  medium. 

6.  id  .aux  is  the  name  of  the  auxiliary  structures. 

7.  sheet.model  is  the  name  of  the  sheet  containing  these  definitions. 

8.  sheet.grid  is  the  name  of  the  sheet  containing  the  definition  of  the  finite 
difference  grid. 

9.  sheet.topo  is  the  sheet  containing  the  definition  of  the  topography. 

10.  sheet.med  is  the  sheet  defining  the  acoustic  medium. 

11.  sheet.aiix  is  the  sheet  defining  the  auxiliary  structures  such  as  buildings, 
paved  areas,  and  bridges. 
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Figure  5.1.  Main  spreadsheet  containing  probiem 
setup. 
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The  id  fields  will  be  strung  together  to  provide  a  unique  name  for  the 
analysis  and  data  files.  This  will  be  called  the  MODELID.  The 
MODELID  is  a  concatenation  of  the  id.model,  id.grd,  id.topo,  and 
id.med  data  items.  For  the  spreadsheet  shown  in  Figure  5.1,  the 
MODELID  is  mctp.  More  than  one  character  can  be  entered  for  an  id 
field,  but  this  would  make  the  name  longer. 

More  than  one  sheet  name  can  be  entered  on  the  line  for  sheet.aux.  For 
example,  if  there  are  eight  sheets  describing  the  buildings,  eight  sheet 
names  can  be  entered  on  this  line  in  the  adjacent  eight  cells.  There  is  an 
option  in  the  auxiliary  sheets  that  instructs  MATLAB  to  modify  the 
topography  to  include  the  tops  of  the  auxiliary  structures.  For  this  to  work 
correctly,  all  auxiliary  structures  must  be  loaded  in  memory  at  one  time. 
Therefore,  all  of  the  sheets  describing  the  structures  must  be  entered  on 
the  sheet.aux  line.  This  allows  the  extraction  of  results  at  a  desired 
distance  above  the  topography  and  the  top  of  the  auxiliary  structures  (e.g., 
the  tops  of  the  roofs  of  buildings).  Otherwise,  all  results  will  be  extracted 
for  a  specific  elevation  or  a  specific  height  above  the  topography.  If  the  flag 
is  used  to  modify  the  topography,  but  each  auxiliary  file  is  run  separately 
by  rerunning  the  MATLAB  processing  routines,  the  results  will  only  be 
output  at  the  top  of  the  auxiliary  structures  for  the  last  sheet  processed. 

The  sheet.topo  can  be  left  blank  to  use  zero  topography,  which  is 
appropriate  for  propagation  without  a  topographic  surface.  An  example 
would  be  acoustic  propagation  well  away  from  a  topographic  surface. 
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All  sheet  names  used  are  case  sensitive.  That  is,  if  a  sheet  is  named 
Buildings  and  buildingS  is  used  on  this  sheet,  an  error  will  occur. 

5.4  Spreadsheet  to  describe  grid 

The  spreadsheet  shown  in  Figure  5.2  contains  information  that  controls 
the  spacing  and  size  of  the  finite  difference  grid.  The  options  are: 

1.  options.grid.type  can  either  be  constant  or  variable.  Constant 

specifies  that  the  spacing  between  the  grid  points  does  not  change  and  is  a 
constant.  Variable  specifies  that  the  spacing  can  vary  by  stretching 
transformations.  Additional  information  is  required  to  define  a  variable 
grid  than  is  shown  in  Figure  5.2. 

2.  grd.dxMin  is  the  spacing  in  the  x  direction. 

3.  grd.dyMin  is  the  spacing  in  the  y  direction. 

4.  grd.dzMin  is  the  spacing  in  the  z  direction. 

5.  grd.nx  is  the  number  of  grid  points  in  the  x  direction. 

6.  grd.ny  is  the  number  of  grid  points  in  the  y  direction. 

7.  grd.nz  is  the  number  of  grid  points  in  the  z  direction. 

A  variable  grid  could  be  used  to  provide  for  fewer  grid  points  in  areas  were 
the  wavelengths  in  the  materials  are  larger  based  on  the  wave  speed  and 
target  frequency  (Ketcham  et  al.  2005).  Since  buildings  covered  the  entire 
area  under  consideration  is  this  analysis,  a  constant  spacing  was  used.  A 
larger  spacing  in  the  z  direction  could  have  been  used  above  the  tops  of  the 
buildings,  but  for  simplicity  a  constant  was  used.  The  data  items 
concerning  a  variable  grid  are  not  shown  in  Figure  5.2. 

Figure  5.2.  Definition  of  finite  difference  grid. 
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5.5  Spreadsheet  to  describe  medium  properties  and  iayers 

The  spreadsheet  shown  in  Figure  5.3  contains  information  that  describes 
the  medium  material  properties  and  air  layering.  The  options  are: 
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1.  op1ions.med.type  sets  the  type  of  material  being  used  for  the  medium. 
The  options  are  elastic,  viscoelastic,  or  poroacoustic.  For  the  acoustic 
analyses,  poroacoustic  was  used. 

2.  options.med.distribution  specifies  how  the  material  properties  are 
described.  The  options  are  layers  or  point. 

3.  med.layer.options.aboveground  specifies  how  the  above-ground 
layers  are  defined.  The  options  are  parallel  and  horizontal. 

4.  med.layer.options.belowground  specifies  how  the  below-ground 
layers  are  defined.  The  options  are  parallel  and  horizontal. 

5.  med.strings.propHeaders  are  the  material  property  headers.  They 
should  remain  as  given. 

6.  med.layer.aboveground.Li  specifies  the  material  properties  of  the 
layer  1  above  the  ground  surface.  The  layers  are  specified  from  the  top 
down. 

7.  med.layer.aboveground.L2  specifies  the  definition  of  layer  2  above  the 
ground  surface. 

8.  med.layer.aboveground.L3  specifies  the  definition  of  layer  3  above  the 
ground  surface. 

9.  med.layer.belowground.Li  specifies  the  definition  of  layer  1  below  the 
ground  surface. 


Figure  5.3.  Definition  of  medium  properties  and  iayers. 
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The  layers  option  defines  the  material  distribution  by  horizontal,  sloped 
(seismic),  or  parallel-to-topography  layering.  The  point  option  defines  the 
material  distribution  using  a  point-cloud  distribution  of  materials,  i.e., 
each  node  in  the  3D  volume  is  assigned  an  integer  value  representing  a 
material. 

The  medium  layers  are  defined  using  either  the  option  parallel  or 
horizontal.  Parallel  will  define  the  layers  parallel  to  the  topography  at 
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the  given  height  while  horizontal  will  define  the  layers  to  be  horizontal  at 
the  given  height. 

The  sheet  shown  in  Figure  5.3  has  three  above-ground  air  layers  and  one 
below-ground  layer.  More  layers  may  be  entered  by  adding  additional  lines 
in  the  sheet  with  a  designation  of  .aboveground.L#  or 
.belowground.L#  where  #  is  the  layer  number.  All  layers  are  defined 
using  a  height  in  meters  above  the  topography. 

The  material  properties  entered  for  the  above  and  below  ground  layers  in 
this  sheet  are  the  default  materials  for  the  finite  difference  mesh.  The 
material  properties  may  be  changed  for  regions  of  the  grid  based  on  the 
material  properties  of  objects  defined  in  other  spreadsheets. 

The  poroacoustic  material  properties  are  explained  in  Wilson  and  Liu 
(2004)  and  given  in  Table  5.1.  The  material  properties  approximate  a 
porous-ground  reflecting/absorbing  surface.  From  Figure  5.3,  the  material 
properties  for  a  poroacoustic  material  for  an  acoustical  analysis  consist  of 
the  density,  wave  speed,  and  optional  flow  resistivity,  porosity,  and 
tortuosity  parameters.  The  air  layers  have  material  properties  consisting  of 
the  density  and  wave  speed.  The  tortuosity  factor  and  heat-capacity  ratio 
properties  were  not  used,  and  are  therefore  blank.  The  tortuosity  factor 
defaults  to  4/3  while  the  heat-capacity  ratio  defaults  to  1.4.  The  density 
and  wave-speed  values  entered  should  be  for  the  fluid  portion  of  the 
porous  material  (i.e.,  the  fluid  that  fills  the  pores  of  the  porous  material). 


Table  5.1.  Porous  material  properties  (Wilson  and  Liu  2004). 


Material 

Flow  Resistivity 

Pa-s/m2 

Porosity 

Tortuosity 

Asphalt 

3.00E+07 

0.1 

3.2 

Grass 

2.00E+05 

0.5 

1.4 

Forest 

l.OOE+05 

0.6 

1.3 

Sand 

5.00E+04 

0.35 

1.6 

Snow 

l.OOE+03 

0.6 

1.7 

5.6  Spreadsheet  to  describe  topography 

The  spreadsheet  shown  in  Figure  5.4  contains  information  that  describes 
the  topography.  The  options  are: 
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1.  options.topo.type  specifies  how  the  topography  is  defined.  The  options 
are  geotiff  and  slope. 

2.  DEM.Filename  is  the  name  of  the  GeoTIFF  file. 

3.  DEM.numZNodesAdjacent  is  the  number  of  nodes  below  the  lowest 
ground  elevation  for  the  acoustic  analysis. 

4.  DEM.iiiitial.gridRotatioii  specifies  a  rotation  for  the  placement  of  the 
GeoTIFF  file.  The  rotation  is  in  degrees  counterclockwise  from  the  positive 
x-axis.  This  value  can  be  changed  in  MATLAB. 

5.  DEM.initial.xEAlign  specifies  the  x  alignment  of  the  GeoTIFF  file.  This 
value  is  the  initial  x  and  Easting  values  that  the  grid  origin  is  aligned  to. 
The  grid  origin  starts  at  x=o.  This  value  can  be  changed  in  MATLAB. 

6.  DEM.iiiitial.yNAligii  specifies  the  y  alignment  of  the  GeoTIFF  file.  This 
value  is  the  initial  y  and  Northing  values  that  the  grid  origin  is  aligned  to. 
The  grid  origin  starts  at  y=o.  This  value  can  be  changed  in  MATLAB. 


Figure  5.4.  Definition  of  topography. 


options.topo.type  ^ 

geotiff 

DEM.Filename 

ElevLIDAR2010-latlon-bridges-2500.tif 

DEM.numZNodesAdjacent  ^ 

20 

DEM. initial. gridRotation  ^ 

DEM.initial.xEAlign  ^ 

DEM. initial. yNAIign  ^ 

The  topography  can  be  defined  by  the  option  slope,  which  would  create 
horizontal  or  sloped  flat  models.  The  geotiff  option  specifies  the 
topography  defined  by  a  GeoTIFF  file.  Global  Mapper  can  output  a 
GeoTIFF  as  described  in  Chapter  3. 

The  DEM.numZNodesAdjacent  option  specifies  the  number  of  nodes 
below  the  lowest  ground  elevation  for  the  acoustic  analysis.  This  number 
of  nodes  should  also  accommodate  the  absorbing  boundary  condition.  The 
thickness  of  the  perfectly  matched  layer  (PML)  is  hardcoded  to  10  nodes  in 
MATLAB.  Therefore,  the  value  of  20  entered  in  the  sheet  includes  10 
nodes  below  the  lowest  elevation  and  10  nodes  for  the  perfectly  matched 
layer.  The  PML  is  the  same  on  all  sides  of  the  mesh. 

The  values  for  the  placement  of  the  GeoTIFF  were  left  blank  and  entered 
when  running  the  MATLAB  preprocessing  routines. 
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5.7  Spreadsheet  to  describe  buildings 

The  spreadsheet  shown  in  Figure  5.5  contains  information  that  describes 

the  buildings.  Multiple  spreadsheets  may  be  used  to  describe  parts  of  the 

total  building  data  set.  The  options  are: 

1.  oplions.aiix.type  specifies  how  auxiliary  structures  are  described.  The 
options  are  vtksDFV,  planForm,  and  landCoverGeoTiff.  For 
buildings,  planForm  was  used. 

2.  aiix.options.xyCoordSys  defines  the  coordinate  system  for  the 
buildings.  This  option  is  either  model  or  UTM.  For  buildings,  UTM  was 
used. 

3.  aiix.options.zCoordSys  defines  how  the  z  elevation  is  defined.  The 
options  are  rel2ModelZero,  rel2SeaLevel,  and  rel2Topo.  For 
buildings,  rel2Topo  was  used. 

4.  aiix.options.modifyTopo  indicates  whether  the  original  topography 
should  be  modified  to  include  the  tops  of  the  buildings.  Yes  or  no  are  the 
options.  For  buildings,  yes  was  used. 

5.  aiix.options.ancliorInterfaces  indicates  whether  the  footprint  of  a 
structure  is  aligned  to  the  FD  grid.  Options  yes  or  no  can  be  chosen.  For 
buildings,  no  was  used. 

6.  aiix.strings.propHeaders  are  the  strings  defining  the  material 
properties  and  should  be  left  as  shown. 

7.  aux.properties.mati  specifies  the  material  properties  of  material  1.  The 
material  entered  represents  asphalt  (Table  5.1). 

8.  a11x.properties.mat2  specifies  the  material  properties  of  material  2.  The 
material  entered  represents  asphalt  (Table  5.1). 

9.  aux.options.type.stmcti  specifies  how  the  structures  are  created.  The 
options  are  rtPrism  and  surfFeat.  For  buildings,  rtPrism  was  used. 

10.  aiix.options.outlineMethod.structi  specifies  the  outline  method  for 
the  structure.  The  options  are  polygon,  rectAsPolygon, 
circleAsPolygon,  or  ellipseAsPolygon.  For  buildings,  polygon  was 
used. 

11.  aux.options.rotationOrigin.structi  defines  the  rotation  point  for  the 
structure.  The  options  are  minBBox  and  centroid.  For  buildings, 
minBBox  was  used. 

12.  aiix.xVertices.stmcti  are  the  x  vertices  of  the  polygon  defining  the 
structure.  The  vertices  are  listed  horizontally  in  cells  adjacent  to  the  data 
option. 
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13.  aiix.yVertices.s1ructi  are  the  x  vertices  of  the  polygon  defining  the 
structure.  The  vertices  are  listed  horizontally  in  cells  adjacent  to  the  data 
option. 

14.  aiix.zHeight.s1ructi  is  the  height  of  the  structure  in  meters. 

15.  aux.zBottom.structi  is  the  bottom  elevation  of  the  structure  in  meters. 

16.  aux.structMatlndex.structi  is  the  number  of  the  material  property 
applied  to  this  structure. 


Figure  5.5.  Definition  of  buiidings. 


options. aux.type 

planForm 

aux.  options.  xyCoordSys 

UTM 

aux.options.zCoordSys  ^ 

rel2Topo 

aux.options.modifyTopo 

yes 

aux.options.anchorlnterfaces  ^ 

no 

aux.repeat.set  ^ 

aux.xyTranslation.set 

aux.rotation.set  ^ 

aux.rotationOrigin.set 

aux.zrotation.set 

aux.  St  rings.  propHeaders 

Material 

Index 

density 

(kg/m3) 

c  (nVs) 

1 

flow 

resistivity 

(Pa-s/m2) 

porosity 

tortuosity 

aux.  properties,  mat  1  ^ 

1 

1.2 

344 

3.00E+07 

0.1 

3.2 

aux.  properties.  nnat2 

2 

1.2 

344 

3.00E+07 

0.1 

3.2 

aux.  options,  type.  St  ructl 

rtPrism 

aux.  options,  out  lineMet  hod.  St  ructl 

polygon 

aux.  options.  rotationOrigin.  St  ructl 

minbbox 

aux.  repeat. structl 

aux.xyT  ranslation. structl 

aux.  rotation.  St  ructl 

aux.xyzSca  ling.  St  ructl 

aux.  wallThickness.  St  ructl 

aux.floorThickness.  St  ructl 

aux.  ceilingThickness.  St  ructl 

aux.xVertices.structl 

708798.4 

708795.8 

708795.7 

708798.3 

aux.yVertices.  St  ructl 

3634920 

3634920 

3634922 

3634922 

aux.  shapeAsPolygon.  St  ructl 

aux.zHeight. structl 

3.3528 

aux.zBottom.structi 

0 

aux.  voidMatIndex.  St  ructl 

aux.structMatlndex.structi 

2 

The  PSTOP3D  program  and  the  MATLAB  processing  routines  refer  to  any 
structures  that  sit  on  or  under  the  topography  as  auxiliary  structures.  The 
option  vtksDFV  allows  3-D  structures  to  be  defined  using  the 
Visualization  Toolkit  (VTK)  format  from  Kitware  (Kitware  2014).  This  is  a 
file  format  for  describing  the  face  and  vertex  data  of  the  3-D  structure.  The 
planForm  option  allows  the  definition  of  a  structure  by  extruding  an 
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outline  in  the  x-y  plane.  Objects  can  either  have  a  thickness  and  z-extents 
relative  to  the  surface  topography  or  extrude  in  the  z-direction  as  right 
prisms.  The  landCoverGeoTiff  option  allows  for  surface  features  to  be 
defined  based  on  a  USGS  land  cover  geoTIFF  file.  The  planForm  option 
is  used  in  this  study  to  model  the  buildings  and  surface  features. 

The  model  option  specifies  coordinates  relative  to  the  model  origin  (x=o, 
y=o),  where  the  UTM  option  specifies  coordinates  in  the  UTM  coordinate 
system.  For  the  acoustic  analyses,  the  UTM  coordinate  system  was  easier 
to  use  since  the  GIS  output  was  in  UTM  coordinates. 

The  rel2ModelZero  option  specifies  the  z  values  relative  to  z=o  (the 
model  origin).  The  rel2SeaLevel  option  specifies  the  z  values  relative  to 
sea  level,  which  is  useful  with  UTM  coordinates  and  known  topographic 
elevations.  The  rel2Topo  option  specifies  the  z  values  relative  to  the 
lowest  topography  value  within  the  object  footprint. 

For  the  acoustic  analyses,  the  rel2Topo  option  worked  best.  Since  a 
height  value  was  specified  for  each  building,  this  height  was  relative  to  the 
topography.  Placing  the  building  at  the  z  elevation  of  the  lowest  elevation 
within  the  building  footprint  ensures  that  the  bottom  of  the  building  is  not 
above  the  terrain  anywhere  within  the  footprint.  This  option  best  models 
the  actual  building  placement. 

The  .modifyTopo  data  item  can  be  used  to  modify  the  original 
topography  to  include  the  tops  of  the  buildings.  PSTOP3D  will  use  this 
modified  topography  when  saving  slice  data  relative  to  the  topography, 
which  is  useful  if  the  results  at  the  tops  of  the  buildings  are  desired.  To 
modify  the  topography  correctly,  all  buildings  must  be  listed  on  the  main 
setup  sheet  as  discussed  in  Section  5.3.  Two  files  are  created  when  this 
spreadsheet  is  processed  by  MATLAB.  The  first  file  ends  with 
_geo.dat.o.origtopo  and  is  the  original  topography  that  has  not  been 
modified,  the  second  file  ends  with  _geo.dat.o.2Doutotpo  and  is  the 
modified  topography. 

The  .anchorlnterfaces  data  item  can  be  used  to  anchor  a  polygon's  x-y 
vertices  at  material  interface  grid  locations;  otherwise,  the  vertices  are 
free-floating.  This  can  be  used  to  make  the  buildings  exactly  line  up  with 
the  FD  grid  nodes.  For  the  acoustic  analyses,  the  vertices  were  left  free- 
floating. 
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The  .type.structi  data  item  for  the  auxiliary  structures  can  be  rtPrism 
(right  prism)  or  surfFeat  (surface  feature).  The  rtPrism  option  indicates 
that  the  structure  should  be  defined  by  extruding  a  polygon  that  defines 
the  footprint  of  the  object.  The  polygon  will  have  a  horizontal  (flat)  bottom 
and  top.  The  surFeat  option  indicates  that  the  object  has  a  designated 
thickness  and  follows  the  topography.  The  buildings  for  the  acoustic 
analyses  used  the  rtPrism  option. 

The  .structi  ending  on  the  data  items  means  that  the  data  item  is  for 
structure  number  i.  Many  structures  can  be  defined  on  one  sheet  and  the 
data  items  would  simply  repeat  themselves  but  with  an  ending  of 
.struct#,  where  #  is  the  number  of  the  structure  being  defined. 

The  .outlineMethod  specifies  how  the  coordinates  of  the  object  are 
defined.  The  polygon  option  allows  the  polygon  vertices  to  be  entered 
directly.  The  other  options  allow  shapes  to  be  used  to  define  the  structure 
coordinates.  Since  the  building  footprint  vertices  are  known,  the  polygon 
option  was  the  easiest  to  use. 

The  .rotationOrigin  data  item  can  be  minBBox  to  use  the  origin  of  the 
bounding  box  for  the  structure  or  centroid  to  use  the  centroid  of  the 
polygon.  For  the  acoustical  analyses,  the  placement  of  the  buildings  was 
from  the  actual  building  layout;  therefore,  buildings  did  not  need  to  be 
rotated.  A  value  is  required,  so  minBBox  was  used. 

The  .zheight  is  the  height  of  the  object  above  the  object's  bottom  for  right 
prisms  or  the  height  of  the  surface  above  the  topography  for  the  surface 
feature.  The  .zheight  will  be  positive  (+)  for  buildings.  The  .zhottom 
data  item  is  used  for  buildings  (right  prisms)  and  is  blank  for  surface 
features. 

5.8  Spreadsheet  to  describe  paved  areas 

The  Urban  Modeler  program  will  output  two  files  for  area  features:  one  file 
will  contain  information  for  the  area  features,  while  the  other  file  will 
contain  information  for  the  islands.  The  spreadsheet  shown  in  Figure  5.6 
contains  information  that  describes  the  paved  areas  representing  roads, 
streets,  highways,  parking  lots,  and  medians.  The  data  items  and  options 
are  similar  to  the  ones  discussed  for  buildings  in  Section  5.7.  The 
differences  in  data  items  are: 
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1.  aiix.op1ions.zCoordSys  is  given  a  value  but  is  not  used  for  paved  areas 
which  are  modeled  as  surface  features. 

2.  aiix.options.modifyTopo  is  set  to  no. 

3.  aux.options.type.structi  is  set  to  surfFeat. 

4.  aux.zHeight.structi  is  set  to  -2  m. 

5.  aux.zBottom.structi  is  left  blank. 


Figure  5.6.  Definition  of  paved  areas. 


options. aux.type 

planForm 

aux.options.xyCoordSys  ^ 

UTM 

aux.options.zCoordSys  ^ 

rel2Topo 

aux.options.  modifyTopo  ^ 

no 

aux.options.anchorlnterfaces  ^ 

no 

aux.repeat.set  ^ 

aux.xyTranslation.set  ^ 

aux.rotation.set  ^ 

aux.rotationOrigin.set  ^ 

aux.zrotation.set 

aux.strings.propHeaders 

Material  Index 

density 
(kg/ m3) 

c  (nVs) 

flow 

resistivity 
(Pa-s/ m2) 

porosity 

tortuosity 

aux.properties.matl  ^ 

1 

1.2 

344 

3.00E-F07 

0.1 

3.2 

aux.properties.mat2 

2 

1.2 

344 

3.00E-F07 

0.1 

3.2 

aux.options.type.structi 

surfFeat 

aux.  options,  out  lineMethod.  St  ructl 

polygon 

aux.options.  rotationOrigin.structl 

minbbox 

aux.  repeat. structl 

aux.xyT  ranslation. structl 

aux.  rotation. structl 

aux.xyzScaling.  structl 

aux.wallThickness.  St  ructl 

aux.floorThickness.  St  ructl 

aux.ceilingThickness.  structl 

aux.xVertices.  structl 

706423.653 

706417.68 

706418.006 

706418.493 

706418.957 

706419.396 

aux.yVertices.  structl 

3636845.629 

3636845.63 

3636847.37 

3636850.11 

3636852.85 

3636855.6 

aux.shapeAsPolygon.  structl 

aux.zHeight.structi 

-2 

aux.zBottom.  St  ructl 

aux.voidMatIndex.  structl 

aux.structMatIndex.  St  ructl 

2 

The  .modifyTopo  data  item  is  set  to  no.  This  item  is  usually  used  for 
buildings. 


The  .type.structi  data  item  for  the  auxiliary  structures  is  set  to  surfFeat 
(surface  feature).  The  surFeat  option  indicates  that  the  object  has  a 
designated  thickness  and  follows  the  topography.  Since  the  paved  areas 
essentially  follow  the  topography,  this  option  adequately  models  the  paved 
features. 


The  .zheight  is  the  height  of  the  surface  above  the  topography  for  a 
surface  feature.  The  .zheight  will  be  positive  (+)  for  surface  features 
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above  the  topography  and  negative  (-)  for  surface  features  below  the 
topography.  The  .zbottom  data  item  is  not  used  for  surface  features  and 
is  left  blank. 

The  GIS  information  defining  the  paved  areas  contains  what  are  known  as 
islands.  The  islands  are  essentially  holes  specified  in  larger  complex 
polygons.  That  is,  a  large  rectangular  polygon  could  be  defined  that 
represents  a  road  network.  Islands  would  model  the  places  where  the 
streets  are  not  present,  essentially  cutting  out  regions  that  would  define  a 
city  block.  Example  islands  are  shown  in  Figure  4.4.  The  process  to  save 
paved  area  information  including  the  islands  is  discussed  in  Section  3.4.5. 

The  definition  of  the  islands  is  the  same  as  the  paved  areas  except  a 
different  material  property  is  used  to  model  the  islands.  Since  the  islands 
are  essentially  a  hole  into  the  underlying  topography,  the  islands  are 
assigned  a  material  property  of  grass  as  given  in  Table  5.1  and  shown  in 
Figure  5.7.  Therefore,  the  data  items  in  Figure  5.7  that  are  different  for  an 
island  versus  a  paved  area  are: 

1.  aiix.properties.mati  specifies  the  material  properties  of  material  1.  The 
material  entered  represents  grass  (Table  5.1). 

2.  a11x.properties.mat2  specifies  the  material  properties  of  material  2.  The 
material  entered  represents  grass  (Table  5.1). 

The  islands  can  be  detected  using  Global  Mapper  as  discussed  in  Section 
3.4.4  and  with  Urban  Modeler  as  discussed  in  Section  4.2.  There  maybe 
differences  between  the  two  sets  of  islands  seen  in  Global  Mapper  versus 
Urban  Modeler.  The  differences  arise  because  of  the  processing  of  the  data 
when  saving  the  paved-area  information  from  Global  Mapper.  If  the  paved 
areas  are  output  using  a  grid  as  discussed  in  Section  3.4.5,  the  features  are 
split  along  the  grid  lines.  Any  islands  that  are  split  by  a  grid  line  will 
disappear.  This  is  because  the  grid  line  separates  the  feature  surrounding 
the  island  into  two  separate  features  that  do  not  require  the  use  of  an 
island. 


ERDC  TR-15-5 


85 


Figure  5.7.  Definition  of  islands  in  paved  areas. 


options. aux.type 

planForm 

aux.options.xyCoordSys  ^ 

UTM 

aux.options.zCoordSys  ^ 

rel2Topo 

aux.options.modifyTopo  ^ 

no 

aux.options.anchorlnterfaces  ^ 

no 

aux.repeat.set 

aux.xyTranslation.set 

aux.rotation.set 

aux.rotationOrigin.set 

aux.zrotation.set 

aux.  strings.  propHeaders 

Material  Index 

density 

(kg/m3) 

c  (nVs) 

flow 

resistivity 

(Pa-s/m2) 

porosity 

tortuosity 

aux.properties.matl  ^ 

1 

1.2 

344 

2.00E+05 

0.5 

1.4 

aux.properties.mat2 

2 

1.2 

344 

2.00E+05 

0.5 

1.4 

aux.  options,  type. struct  1 

surfFeat 

aux.  options.  outlineMethod.  struct  1 

polygon 

aux.  options.  rotationOrigin.  struct  1 

minbbox 

aux.  repeat. structl 

aux.xyT  ranslation. structl 

aux.  rotation. structl 

aux.  xyzScaling.  structl 

aux.wallThickness.  structl 

aux.floorThickness.  structl 

aux.  ceilingThickness.  structl 

aux.  xVertices.  structl 

706555.167 

706401.398 

706396.036 

706394.636 

706393.257 

706391.899 

aux.  yVert  ices,  structl 

3637095.697 

3637090.86 

3637088.72 

3637086.47 

3637084.21 

3637081.94 

aux.  shapeAsPolygon.  structl 

aux.zHeight. structl 

-2 

aux.zBottom.structi 

aux.  voidMatIndex.  structl 

aux.  structMatIndex. structl 

2 

5.9  Spreadsheet  to  describe  bridges 

The  spreadsheet  shown  in  Figure  5.8  contains  information  that  describes 
the  bridge  objects.  The  data  items  and  options  are  similar  to  the  ones 
discussed  for  buildings  in  Section  5.7.  The  differences  in  data  items  are: 


1.  aux.options.zCoordSys  is  set  to  rel2SeaLevel. 

2.  aiix.zHeight.structi  is  the  height  of  the  structure  in  meters. 

3.  aux.zBottom.structi  is  the  bottom  elevation  of  the  structure  in  meters. 


The  rel2SeaLevel  option  specifies  the  z  values  relative  to  sea  level,  which 
is  useful  with  UTM  coordinates  and  known  topographic  elevations.  For 
bridges,  the  Urban  Modeler  program  computes  a  minimum  elevation  of 
the  bridge  based  on  the  z-coordinates  assigned  to  the  bridge  vertices  by 
Global  Mapper  as  discussed  in  Section  3.5.1.  This  minimum  elevation  is 
the  .zbottom  data  item.  Since  this  elevation  is  in  UTM  coordinates,  the 
rel2SeaLevel  option  works  best. 
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Bridges  use  the  rtPrism  (right  prism)  option  as  described  Section  5.7. 
The  .zheight  is  the  height  of  the  object  above  the  object's  bottom  for  right 
prisms.  For  bridges,  .zheight  is  assigned  a  value  of  2  m. 


Figure  5.8.  Definition  of  bridge  objects. 
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planForm 

aux.  opt  ions.  xyCoordSys 
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aux. options. zCoordSys  ^ 
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aux.  opt  ions.  type.  St  ructl 

rtPrism 

aux.options.outlineMethod.structl 

polygon 

aux.options.rotationOrigin.  St  ructl 

minbbox 

aux. repeat. structl 

aux.xyT  ranslation. structl 

aux. rotation. structl 

aux.xyzSca  ling.  St  ructl 

aux.wallThickness.  St  ructl 

aux.floorThickness.  St  ructl 

aux.ceilingThickness.  St  ructl 

aux.xVert  ices.  St  ructl 

707910.163 

707910.205 

707913.072 

707915.297 

707918.473 

707922.016 

a  ux.y  Vert  ices.  St  ructl 

3635397.488 

3635397.42 

3635394.68 

3635393.42 

3635392.38 

3635392.13 

aux.shapeAsPolygon.structl 

aux.zHeight. structl 

2 

aux.zBottom.structl 

171.9 

aux.voidMatlndex.structl 

a  ux.  struct  Matindex.  St  ructl 

2 

aux.  repeat. structl 
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6  Run  MATLAB  to  Produce  PST0P3D  Data 
Files 

6.1  General 

Once  the  Excel  spreadsheets  are  composed,  the  MATLAB  preprocessing 
routines  can  be  run.  This  corresponds  to  Step  6  on  the  process  diagram  in 
Figure  1.4.  To  run  the  MATLAB  routines: 

1.  Execute  MATLAB  in  Windows. 

2.  On  the  main  MATLAB  screen  shown  in  Figure  6.1,  navigate  to  the 
directory  where  the  scripts  are  installed.  The  main  script  to  execute  is 
called  rvgModel.m. 

3.  The  path  may  need  to  be  changed  to  include  the  directory  and  sub¬ 
directories  where  the  scripts  are  installed.  On  the  Home  tab,  click  the  Set 

Path  button. 

4.  As  shown  in  Figure  6.2,  click  Add  with  Subfolders  and  select  the 
appropriate  directory.  Click  Save. 

5.  Double-click  on  the  rvgModel.m  script  shown  in  Figure  6.1. 

6.  The  editor  with  rvgModel.m  will  appear  as  shown  in  Figure  6.3.  Click  on 
the  Run  button  on  the  main  toolbar. 

7.  The  main  menu  for  the  processing  routines  will  appear  as  shown  in 
Figure  6.4. 

8.  On  the  menu,  select  the  button  Select  spreadsheet  with  input  data. 

9.  Navigate  to  where  the  spreadsheet  is  saved  and  click  on  the  file. 

The  menu  in  Figure  6.4  references  the  Excel  1995  format.  If  the 
spreadsheet  is  in  the  Excel  1995  format,  the  spreadsheet  can  only  contain 
16,384  rows  (lines).  The  amount  of  data  that  needs  to  be  stored  in  the 
spreadsheets  exceeded  these  limits.  Since  the  processing  was  performed 
on  a  Windows  7  machine,  the  Excel  2007  format  was  a  better  format  to 
use.  The  Excel  2007  format  allows  1  million  lines  per  spreadsheet. 

The  processing  of  the  spreadsheets  is  performed  in  the  order  of  the 
buttons  on  the  menu  in  Figure  6.4.  Therefore,  click  the  buttons  and 
perform  the  processing  from  the  top  of  the  menu  downward. 


Ift 
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Figure  6.1.  Main  MATLAB  screen. 


Figure  6.2.  Set  path  for  MATLAB. 
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Figure  6.3.  Running  rvgModel.m. 


Figure  6.4.  Main  menu  for  MATLAB  preprocessing  routines. 


y  Main  menu 


J2 


-Select  step 


Clear  existing  figures 


Move  to  folder  with  output  files 


Select  spreadsheet  with  input  data 


Generate  grid 


Generate  topography  function 


Define  propagation  medium 


This  code  supports  the  parallel 
FDTD  code  pstopSd.  It  creates  a 
model  for  seismic  or  acoustic 
propagation  calculations. 

Input  is  primarily  through  Excel 
spreadsheets  in  Excel  95  format,  or 
files  from  other  codes  (such  as 
Open  Office)  saved  as  such. 

S.  Ketcham,  USAERDC, 

Stephen.A.Ketcham@usace.army.mi 

I 


Define  auxiliary  structure 


Generate  source 


Finish 


ERDC  TR-15-5 


90 


6.2  Generate  grid 

To  generate  the  grid,  click  the  Generate  grid  button  shown  in  Figure  6.4. 
The  plots  shown  in  Figures  6.5  through  6.7  will  display.  These  plots 
display  the  variable  grid  parameters  for  the  grid,  which  is  the  relationship 
between  the  stretching  transformations  and  the  physical  grid.  The 
relationship  is  seen  to  be  linear  from  the  figures  and  the  derivative  is  a 
constant  value  of  1,  denoting  a  one-to-one  relationship  between  the 
transformed  grid  and  the  physical  grid.  These  parameters  describe 
constant  grid  spacing. 

After  the  grid  is  processed,  the  informational  dialog  shown  in  Figure  6.8  is 
displayed  that  reports  the  limits  of  the  grid. 

Figure  6.5.  Variable  grid  parameters  denoting  a  constant  grid  in  the  x-direction. 
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Figure  6.6.  Variable  grid  parameters  denoting  a  constant  grid  in  the  y-direction. 
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Figure  6.7.  Variable  grid  parameters  denoting  a  constant  grid  in  the  z-direction. 
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Figure  6.8.  Informational  dialog  displaying  grid 
limits. 


BB  Current  Information  [' 

r 

Length  of  model  in  x  direction  =  nx  =  249G. 

Length  of  model  in  y  direction  =  ny  =  249G. 

Length  of  model  in  z  direction  =  nz  =  25G. 

The  model  output,  without  translation  and  rotation,  will  span: 

X  =  0  to  2498.9995  m,  and 

V  =  0  to  2498.9995  m,  and 
z  =  0  to  255  m. 

1 

k. 

6.3  Generate  topography  function 

To  generate  the  topography: 

1.  Click  the  (Generate  topography  function  button  shown  in  Figure  6.4. 
The  plots  shown  in  Figure  6.9  along  with  the  informational  dialog  shown 
in  Figure  6.10  will  display. 

2.  Figure  6.11  will  display  and  the  x  (Easting)  and  y  (Northing)  location  of  the 
grid  can  be  input  as  shown. 

3.  Click  on  the  Draw  grid  on  DEM  button  and  Figure  6.12  will  display 
showing  the  location  of  the  grid  on  the  topography.  The  dialog  in 
Figure  6.13  will  display  with  information  about  the  location  of  the  grid. 


Figure  6.9.  Display  of  topography  as  defined  in  spreadsheet. 
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Figure  6.10.  Informational  dialog  giving 
information  about  topography. 


Figure  6.11.  Defining  the 
position  of  the  grid  on  the 
topography. 


Figure  6.12.  Display  of  location  of  grid  on  topography. 
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Figure  6.13.  Informational  dialog  about 
grid  location. 


6.4  Define  propagation  medium 

To  define  the  propagation  medium: 

1.  Click  on  the  Define  propagation  medium  button  in  Figure  6.4. 

2.  The  dialog  to  generate  slice  plots  will  display  as  shown  in  Figure  6.14. 
Choose  a  slice  location  by  moving  the  sliders  and  click  Draw  figure(s). 

3.  The  slice  plot  will  display  as  shown  in  Figure  6.15.  The  plot  may  be  zoomed 
in  and  panned  to  display  the  data  as  desired.  This  plot  shows  four  layers, 
which  corresponds  to  the  data  shown  in  Figure  5.3.  The  green  line  in 
Figure  6.15  denotes  the  top  of  the  topography. 

4.  Click  Finish. 


Figure  6.14.  Input  for  slice  plots. 


ERDC  TR-15-5 


95 


Figure  6.15.  Slice  showing  the  layers  in  the  medium. 


6.5  Define  auxiliary  structures 

To  define  the  auxiliary  structures: 

1.  Click  on  the  Define  auxiliary  structure  button  in  Figure  6.4. 

2.  The  spreadsheets  listed  on  the  sheet.aux  line  of  the  main  spreadsheet 
discussed  in  Section  5.3  will  be  processed.  A  single  sheet  or  multiple  sheets 
maybe  processed.  If  a  single  sheet  is  processed,  the  name  of  the 
spreadsheet  on  the  sheet.aux  line  can  be  changed,  the  spreadsheet  saved, 
and  the  MATLAB  rvgModel.m  script  run  again.  All  previously  processed 
data  will  be  read  in  and  the  additional  auxiliaiy  sheet  will  be  processed. 
This  can  be  done  as  many  times  as  required  to  process  the  auxiliary 
structures  sheet  by  sheet. 

3.  Figures  6.16  and  6.17  show  the  result  of  processing  partial  building  data. 

4.  After  a  sheet  containing  auxiliary  structures  is  processed,  the  slice  plot 
selection  dialog  as  shown  in  Figure  6.18  is  displayed. 

5.  Move  the  sliders  to  select  a  slice  location.  Click  Draw  figure(s). 

6.  Slice  plots  are  displayed  for  the  selected  locations.  The  x-z  plane  slice  plot 
is  shown  in  Figure  6.19.  The  green  line  in  the  figure  represents  the 
topography.  The  red  squares  represent  nodes  defining  an  auxiliaiy 
structure.  As  seen  from  Figure  6.19,  the  topography  has  been  modified  to 
include  the  tops  of  the  buildings. 

7.  Figures  6.20  through  6.23  show  various  sets  of  processed  buildings. 
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8.  Figure  6.24  shows  the  processed  paved  areas  for  a  strip  that  represents 
one-fifth  of  the  paved  areas.  This  set  does  not  include  islands  and  is 
defined  with  a  material  property  of  asphalt  as  discussed  in  Section  5.8. 

9.  Figure  6.25  shows  the  processed  paved  areas  for  the  islands  that  exist  for 
the  areas  defined  in  Figure  6.24.  These  areas  are  defined  with  a  material 
property  of  grass. 

10.  Click  Finish  to  exit  the  dialog  shown  in  Figure  6.18. 


Figure  6.16.  Elevation  data  for  partial  building  data. 


Figure  6.17.  Partial  building  data  shown  in  3-D. 
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Figure  6.18.  Slice  plot  selection  for  auxiliary  structures. 


Figure  6.19.  x-z  plane  slice  plot  showing  buildings. 
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Figure  6.20. 3-D  plot  of  processed  aggregated  buildings  close  to  U.S.  Highway  75. 
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Figure  6.21.  3-D  plot  of  processed  buildings  greater  than  51  ft  in  height. 


Figure  6.22.  3-D  view  of  processed  individual  buildings. 
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Figure  6.23.  Magnified  3-D  view 
of  processed  buiidings. 
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Figure  6.24.  Processed  partial  paved  areas 
without  islands. 
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Figure  6.25.  Processed  partial  islands  of  paved 
areas. 
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6.6  Generate  source 

The  MATLAB  script  rvloadset_io.m  is  run  to  define  a  source  and  its  type, 
location,  and  magnitude.  The  script  is  run  similar  to  the  rvgModel.m  script 
described  in  Section  6.1.  It  also  creates  several  files:  a  header  file  needed  to 
compile  PSTOP3D,  a  file  defining  the  run  parameters  for  PSTOP3D,  and 
two  files  that  define  the  source.  These  files  are  discussed  in  more  detail  in 
Chapter  7. 

The  question-and-answer  sequence  is  shown  in  the  following  text  in 
Courier  font.  The  user  may  press  the  Enter  key  to  accept  the  default  text 
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located  between  the  left  and  right  chevrons  (<  and  >).  Green  text  denotes 
an  input  response.  Red  text  is  additional  explanation. 


C : \rvloadset_20120606\rvloadset_10 
16-NOV-2012;  computer  type:  PCWIN64 

-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k 

POC :  S.  Ketcham,  USAERDC,  Stephen.A.Ketcham@erdc.usace.army.mil. 

rvloadset_V.m  generates  a  set  of  inputs  to  be  applied  to  a  geo_model  in  PST0P3D. 
For  each  input  the  user  sets: 

(1)  Input  location 

(2)  Input  amplitude 

(3)  Input  time  series  (the  signal  as  a  function  of  time 

For  seismic  simulations,  each  input  is  a  set  of  force  vectors  converted 
to  body  forces  using  the  cell  volume  in  the  3D-d  grid. 

For  acoustic  simulations,  each  input  is  an  acoustic  volumetric  source 
at  the  input  location. 

rvloadset_V.m  uses  the  "geo  model"  file  generated  by  rvgeo_V.m. 

rvloadset_V.m  generates  PST0P3D  input  files  COMMAND. FIL  and  pstop3d.h.  Open  these 
files  in  a  text  editor  &  read  instructions  for  preparing  them  for  final  use. 

The  PST0P3D  m-files  should  be  in  your  matlabpath.,  and  Matlab's  "current 
directory" 

is  where  rvgeo_V.m  "geo  model"  will  be  read,  and  where  rvloadset_V. m  output  will 
be  written. 

Clear  all  figures  (y/n) ?  <  n  >  :  y 
INPUT  "GEO"  MODEL  FILES 

-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k 


Input  ID  of  previously  generated  "geo"  model 
(i.e.,  the  characters  before  "_geo.dat")  <  >  :  fctj 

This  name  is  generated  based  on  the  information  given  in  the  spreadsheet  in  Figure 
5.1  and  discussed  in  Section  5.3.  This  ID  is  the  MODELID  and  is  a  concatenation  of 
the  id. model,  id. grid,  id.topo,  and  id.med  data  items. 


rvgeo_readgeo . m : 


Binary  data  from  "geo"  model  files. 


Beginning  reading  "geo"  model  data  from  mctp_geo.dat  files. 

Finished  reading  grid  and  model  parameters  from  mctp_geo . dat . 4 . 

Finished  reading  DEM  origin  parameters  from  mctp_geo . dat . 5 . 

Finished  reading  x,y,z  grid  vectors  from  mctp_geo . dat . 3 . 

Finished  reading  topography  elevation  function  from  mctp_geo . dat . 0 . 

Finished  reading  topography  elevation  function  from  mctp_geo . dat . 0 . origtopo, 
which  is  used  by  PST0P3D  to  save  3D  output  relative  to  topography  without 
auxiliary  structures. 

Finished  reading  indexed  material  properties  from  mctp_geo . dat . 1 . 

Finished  reading  layer  elevation  functions  from  mctp_geo . dat . 01  (zO 

Finished  reading  layer  elevation  functions  from  mctp_geo . dat . 02  (zO 

Finished  reading  layer  elevation  functions 

Finished  reading  layer  elevation  functions 

Finished  reading  layer  elevation  functions  from  mctp_geo.dat 

Completed  data  read  from  "geo"  files  mctp_geo.dat. 


from  mctp_geo.dat. 
from  mctp_geo.dat. 


03 

04 

05 


(zO 
(zO 
( zO 


for  layer 
for  layer 
for  layer 
for  layer 


for  layer  5) 


Size  of  2D  topography  matrices  "z0_topo"  returned  by  rvgeo_readgeo . m  is: 
NX  MODEL  SIZE+2*LO  x  NY  MODEL  SIZE+2*LO  =  2498  x  2498. 


Size  of  all-layers  topography  matrix  "z0_layers"  returned  by  rvgeo_readgeo . m  is: 
num_layers  +  1  x  NX_MODEL_SIZE+2*LO  x  NY_MODEL_SIZE+2*LO  =5x2498x  2498. 

For  example:  z0_layers  includes  the  upper  grid  boundary  max(z)  in 
z0_layers ( 1 , : , : ) ;  surfaces  of  increasingly  lower  layers  follow  in 
z0_layers (2, : , : ) ,  z0_layers (3, : , : ) ,  etc;  z0_layers (num_layers+l , : , : ) 
includes  the  lower  grid  boundary  min(z) . 


3D  material  indices  matrix  "matl  indices"  not  returned;  obsolete. 
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Length  of 
Length  of 
Length  of 
Length  of 
Length  of 
Length  of 
Number  of 


X  array  returned  by  rvgeo_readgeo . m  is  NX_M0DEL_SIZE+2*L0 . 
y  array  returned  by  rvgeo_readgeo . m  is  NY_M0DEL_SIZE+2*L0 . 
z  array  returned  by  rvgeo_readgeo . m  is  NZ_M0DEL_SIZE+2*L0 . 
x_stag  array  returned  by  rvgeo_readgeo . m  is  NX_M0DEL_SIZE+2*L0-1 . 
y_stag  array  returned  by  rvgeo_readgeo . m  is  NY_M0DEL_SIZE+2*L0-1 . 
z_stag  array  returned  by  rvgeo_readgeo . m  is  NZ_M0DEL_SIZE+2*L0-1 . 
materials  returned  by  rvgeo_readgeo . m  is  4. 


rvgeo_readgeo . m  complete. 


ASSIGN  ANALYSIS  ID  AND  PARALLEL-TO-TOPOGRAPHY  OUTPUT  LAYER 

-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k 


Enter  a  short  descriptive  analysis  ID  to  append  to  mctp 
(suggestion:  use  <=  4  characters  to  distinctly  describe  source)  <  >  : 

This  descriptor  for  the  analysis  is  the  LOADID.  This  ID  and  the  previously  entered 
"geo"  ID  will  be  used  when  naming  the  load  files  created. 

Enter  number  of  nodes  above  topo  surface  of  output  layer  that  will  be 
draped  over  topo  surface  (can  be  changed  later  by  editing  COMMAND . FIL) ; 
for  example,  use  -1  to  get  seismic  surface  response,  or  10  to  get 
acoustic  response  at  approximately  one  shortest  wavelength  above  ground 

<10>  :  H 

Using  1  node  will  extract  results  1  node  above  the  topography.  If  the  modifyTopo 
option  was  used  as  described  in  Section  5.7,  then  the  results  will  be  extracted  1 
node  above  the  topography  and  1  node  above  the  tops  of  the  buildings. 

The  lowest  DT  from  the  Courant  condition,  found  in  mctp_geo.txt 
or  auxiliary-file  text  outputs  for  mctp,  is  the  estimated 
governing  criterion  for  time-step  stability.  Enter  the  DT  to  be  used 
for  the  pstop  analysis  (suggestion,  use  0.995  x  the  lowest  DT) 

<  >  :  10.00051 

This  value  must  be  selected  to  satisfy  accurate  representation  of  the  wavelength 
of  interest  and  also  to  provide  stability  of  the  solution.  This  is  discussed  in 
Chapter  7 . 

LOAD  NUMBER  1:  INPUT  BASIC  SOURCE  TYPE 

-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k 


Select  from  the  following  basic  source  types,  which  will  be  given  a  time  series 

history 

1.  "BODY  FORCE,"  force  vector  components  with  selected  configuration 
and  distribution  around  source  location  point. 

2.  "CONSTRAINED  PARTICLE  VELOCITY,"  particle  vector  vector  components 
with  selected  configuration  and  distribution. 

3.  "Q  SOURCE,"  volume  source,  i.e.,  expressed  as  a  dilatation  rate, 
at  a  source  location  point. 

Note:  additional  functions  can  be  added  to  the  program. 

Input  number  of  basic  source  type,  to  be  refined  with  later  inputs:  : 

LOAD  NUMBER  1:  INPUT  Q- SOURCE  DATA 

-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k 


LOAD  NUMBER  1 :  INPUT  SOURCE  LOCATION 

-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k 

»»»»»  rvloadset_location  .m:  Specification  of  source-center  location. 
fnam_command_aux_loc_prompt  assigned  as  global  variable  in  rvloadset_location . m 
Input  loads  graphically?  (y/n)  <  n  >  :  | 

Default  source  location  is  at  mid-domain  (m) : 

Input  source  location  [x,y]  (m)  <  1250.0005  1250.0005>  :  |2 4 7 0  2 iTo] 
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This  is  the  location  of  the  source.  Since  the  analysis  region  is  2500  meters 
square,  the  source  is  located  30  meters  in  from  the  upper  right  corner  of  the 
region . 

Source  1  x-location  =  2469.9529  m  at  x-dir  node  #  2468 

Source  1  y-location  =  2469.9529  m  at  y-dir  node  #  2468 

At  this  location  the  land  elevation  =  49.7091  m  at  z-dir  node  #  52 

At  this  location  the  output  layer  =  51  m  at  z-dir  node  #  53 

Input  desired  source  height  (default  =  output  layer)  (m)  j5^|  : 

Source  1  z-location  =  51  m  at  z-dir  node  #  53 

Do  you  want  to  keep  this  source  location?  (y/n)  [in 

The  center  of  the  source  is  x=2469.9529,  y=2469.9529,  z=51  m. 

<<<<<<<<<<  rvloadset_location . m  complete. 

>>>>>>>>>>  rvloadset_Qsourcedata .m:  Q-source  magnitude  specification  and  output  to 
file . 


The  source  term  in  the  PDE  describing  acoustic  propagation  is  the  applied  change 
in 

volume  per  unit  time  per  unit  volume,  given  the  symbol  "Q"  and  has  units  m^3/s/m^3 
=  1/s. 

Since  we  usually  think  in  terms  of  sound  pressure  (Pa) ,  we  present  three  choices 
to  go 

from  pressure  to  dilatation-rate  time  series  acting  on  a  single  3D  grid  element: 

The  final  source  history  will  be  Q*g(t),  where  g(t),  selected  later,  is  a  unitless 
time  series. 


1.  "NOMINAL  PRESSURE"  calculates  Q(t)  =  (nominal  pressure) / (bulk  modulus^dt ) *g ( t ) 
where  the  nominal  pressure  is  a  realistic  value  for  linear  computations. 

2.  "MEASURED  PRESSURE,"  calculates  Q(t)  =  (4*pi/density) *int [ (P*g (t) @lm*dt] , 
where  the  input  pressure  P(t)  is  assumed  to  be  measured  at  1  m  from  the  source 
in  free  space,  and  numerical  integration  is  done  to  calculate  Q(t) 

3.  "ANALYTICAL  PRESSURE,"  calculates  Q(t)  =  (4*pi/density) *int [ (P*g (t) @lm*dt] , 
where 

the  input  pressure  P(t)  is  an  analytical  function,  and  Q(t)  includes 
the  analytical  integration  of  this  function  (e.g.,  FHarm_Pulse,  FInt_Harm,  etc.) 
Note:  additional  source  functions  and  source  types  can  be  added  to  this  program. 


Input  number  of  selected  Q  source  type: 


The  source  will  be  defined  using  a  nominal  pressure. 

Enter  nominal  pressure  amplitude  (Pa)  to  scale  time  series  |C1000>|  : 

For  the  4  acoustic  materials  of  the  "geo"  model,  the  bulk  moduli  are: 

1:  c  =  331  m/s;  density  =  1.13  kg/m3;  bulk  modulus  =  123803.9295  Pa. 

2:  c  =  337  m/s;  density  =  1.18  kg/m3;  bulk  modulus  =  134011.414  Pa. 

3:  c  =  344  m/s;  density  =  1.2  kg/m3;  bulk  modulus  =  142003.2056  Pa. 

4:  c  =  179.8444  m/s;  density  =  6.272  kg/m3;  bulk  modulus  =  202861.6969  Pa. 


Enter  the  material  bulk  modulus  at  the  source  point  (consider  auxiliary  material 
properties  if  source  is  within  an  auxiliary  material;  they  are  not  listed  above) 
<123803. 9295>  :  |142003.20'^ 

For  this  nominal  pressure  and  bulk  modulus,  dilatation-rate  amplitude  Q  =  14.0842 
(1/s)  . 


Writing  acoustic  source  amplitude  &  location  to  PSTOP3D  file  mctp_UR15_ls.dat  ... 
This  is  one  of  the  files  needed  for  the  source  definition. 


<<<<<<<<<<  rvloadset_Qsourcedata .m  complete. 

LOAD  NUMBER  1:  INPUT  SOURCE  HISTORY 

************************************************************************** 


>>>>>>>>  rvloadset_history .m:  Information  &  specification  of  source  time-series. 

history_mf ilename_prompt  f ilter_f lag_prompt  delay_f lag_prompt  delay_prompt 
nomalize_gt_f lag_prompt 


ERDC  TR-15-5 


103 


assigned  as  global  variable  in  rvloadset_history . m 

Acoustic  wave  speeds  and  1-km  travel  times: 
c  (m/s)  = 

[331.0000  337.0000  344.0000] 

Duration  (s)  for  acoustic  wave  to  travel  1000  m  = 

[3.0211  2.9674  2.9070] 

Estimated  accurate  analysis  bandwidth  using  10  nodes  per  minimum  acoustic 
wavelength  as  the  criterion  to  mitigate  grid  dispersion  errors: 

For  acoustic  wavelength  =  lO^max ( [Dxi, Dkappa, Deta] )  =  10.016  m, 
the  maximum  allowable  frequency  of  the  slowest-speed  material  is 
33.047  Hz  (this  is  a  reasonable  governing  bandwidth) . 

The  maximum  allowable  frequencies  of  the  next-slowest-speed  materials  are 
33.6461  34.3449  Hz 

Time-series  generating  functions  (rvloadset_window_* . m)  in  PSTOP3D  m-files 
directory 

C : \Inf rasound\Current  Version\rvloadset_20120606 : 
rvloadset_window_FC4pulse . m 
rvloadset_window_FGauspuls .m 
rvloadset_window_FHarm_Pulse . m 
rvloadset_window_FHarmonic . m 
rvloadset_window_FInteg_Harm . m 
rvloadset_window_FPulse . m 
rvloadset_window_FRicker . m 
rvloadset_window_FexGauspuls .m 

Input  root  filename  of  time-series  "window"  function  that  defines  source  history 
(this  window  acts  on  previously  input  source  to  give  full  definition  of  source 
magnitude  vs  time--the  function  must  be  written  with  correct  i/o  format) 

<  >  I  rvloadset  window  FRick^ 

This  will  define  a  Ricker  pulse.  Other  options  are  a  Gaussian  pulse  and  a  harmonic 
pulse . 


>>>>>>>>>>  rvloadset_window_FRicker .m:  Ricker-pulse  source  window  function. 

f_ctr_prompt  and  tp_prompt  assigned  as  global  variables  in 
rvloadset_window_FRicker . m 

Enter  center  frequency  (Hz)  for  Ricker  pulse  [50] : 

This  is  the  frequency  of  the  source  that  is  desired. 

Enter  greater  time  (s)  to  move  Ricker  pulse  peak  to  later  time 
(default  starts  Ricker  pulse  at  0  s)  [[0 . 07  451]: 

<<<<<<<<<<  rvloadset_window_FRicker . m  complete. 

>>>>>>>>>>  rvloadset_gtplot .m:  Source  time-series  plot  with  Fourier  and  energy 
spectra . 


<<<<<<<<<<  rvloadset_gtplot . m  complete. 

From  Figure  1, 

at  f  =  15  25  28  33  41  48  Hz,  the  source  energy 

is  down  0  6  10  20  40  60  dB,  respectively. 

The  plots  of  the  Ricker  pulse  are  displayed  in  Figure  6.26.  Figure  6.26a  is  the 
unitless  time  variation  of  the  Ricker  pulse,  Figure  6.26b  is  the  normalized 
Fourier  amplitude  of  the  pulse,  and  Figure  6.26c  is  the  energy  of  the  pulse  in  dB 
relative  to  the  peak  energy.  The  numbers  above  represent  the  circles  on  the  graphs 
in  Figure  6.26.  These  numbers  should  be  used  to  compare  the  source  bandwidth  and 
the  analysis  bandwidth  (that  can  be  accurately  represented)  with  the  bandwidth 
associated  with  the  10  node  per  wavelength  criterion.  The  decibels  for  the  source 
energy  should  be  -40  dB  or  below  for  the  maximum  frequency  that  can  be  represented 
or  a  low-pass  filter  will  need  to  be  applied  to  the  source. 
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Do  you  want  to  keep  this  signal?  (y/n)  [y]  | 

IMPORTANT:  Confirm  source-  signal-  spectral  amplitude  is  not  significant 
above  accurate-analysis-bandwidth  frequency.  Judge  based  on  estimated 
accurate-analysis  bandwidth  above. 

Low-pass  filtering  the  time-  series  signal  g(t)  lowers  high-frequency  amplitude. 
Do  you  want  to  lowpass  filter  g(t)?  (y/n)  <  y  >  :  | 

Adding  a  delay  to  the  time-  series  signal  g(t)  starts  the  signal  at  a  later  time. 

Do  you  want  to  add  a  delay?  (y/n)  |<  n  >|  : 

Normalizing  the  time-  series  signal  g(t)  make  the  maximum  absolute  value  =  1, 
in  which  case  the  source  amplitude  will  be  the  previously  entered  value. 

Do  you  want  to  normalize?  (y/n)  Jn^  : 

>>>>>>>>>>  rvloadset_gtplot .m:  Source  time-series  plot  with  Fourier  and  energy 
spectra . 


<<<<<<<<<<  rvloadset_gtplot . m  complete. 

++++++++++  rv_print.m 

mctp_UR15_loadl_gt . f ig  written  to  disk. 
mctp_UR15_loadl_gt.jpg  written  to  disk. 

++++++++++  rv_print.m  complete. 

Writing  g(t)  to  PST0P3D  file  mctp_UR15_gt.dat 

This  is  the  second  file  needed  for  the  definition  of  the  source.  The  g(t)  function 
is  a  normalized  unit-less  variation  of  the  source  with  time. 


rvloadset_history . m  complete. 


LOAD  NUMBER  1:  OUTPUT  TO  C0MMAND.FIL.mctp_UR15.txt  (load  #1  only) 
and  PARAMETER  FILE  pstopSd . h . mctp_URl5 . txt 

-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k 

»»»»»  rvloadset_commandfile  .m:  Output  of  "COMMAND"  file  and  parameter  file. 
Writing  C0MMAND.FIL.mctp_UR15.txt  ... 

This  file  contains  information  to  run  PSTOPSD.  The  name  is  based  on  the  previous 
input  for  the  "geo"  model  name  and  the  4  character  analysis  description. 

Writing  pstopSd . h . mctp_UR15 . txt  ... 

This  file  contains  the  header  information  needed  to  compile  PSTOPSD.  The  name  is 
based  on  the  previous  input  for  the  "geo"  model  name  and  the  four-character 
analysis  description. 

<<<<<<<<<<  rvloadset_commandf ile . m  complete. 

Add  another  load  to  the  load  set?  (y/n)  <  n  >  :  | 

EINAL  OUTPUT  TO  pstopSd . h . mctp_UR15 . txt 

icicicicicicicicicicicicicicicicicicicicicicicicicicicicicicicicicicic-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k 


rvloadset_10 . m  complete. 
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Figure  6.26.  Definition  of  a  Ricker  pulse. 


Source  time  series  Source  frequency  content 


a)  Time  variation  of  pulse  b)  Fourier  amplitude  of  pulse 


Source  energy  spectrum 


c)  Pulse  energy  in  dB  relative  to  peak  energy 
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7  Compile  and  Run  PSTOP3D 

7.1  PST0P3D  discussion 

This  chapter  discusses  specifics  about  the  PSTOP3D  program  including 
files  necessary  to  compile  and  run  the  program.  This  corresponds  to  Step  7 
of  the  process  diagram  shown  in  Figure  1.4. 

7.1.1  General  background 

The  theoretical  background  for  the  PSTOP3D  code  can  be  found  in  Wilson 
and  Liu  (2004)  and  Ketcham  et  al.  (2005).  Some  of  the  capabilities  of 
PSTOP3D  are  (Ketcham  2006): 

1.  PSTOP3D  allows  the  acoustic  propagation  in  static  medium,  i.e.,  the 
acoustic  propagation  is  not  coupled  with  wind. 

2.  PSTOP3D  uses  a  rectangular  variable  staggered  grid.  The  spacing  of  the 
grid  can  be  a  constant  value  or  may  vary  based  on  stretching 
transformations.  The  spacing  of  the  physical  grid  is  variable  while  the 
spacing  of  the  computational  grid  is  a  constant.  The  staggered  grid  results 
in  more  efficient  computations. 

3.  PSTOP3D  is  a  parallel  processing  code  that  uses  the  message  passing 
interface  (MPI)  software  library  for  parallelization  by  spatial-domain 
decomposition  and  message  passing  at  the  edges  of  the  domain. 

4.  The  finite  difference  time-domain  implementation  uses  second-order 
finite  differences,  which  spatially  allow  highly  discontinuous  material 
interfaces  and  temporally  allows  time  steps  just  below  the  Courant  time 
step. 

5.  The  second-order  acoustic  implementation  has  the  ability  to  model  ground 
surface  and  structures  in  the  model  by  directly  modeling  their  contact  with 
air.  The  code  uses  a  porous-material  model  discussed  in  Wilson  and  Liu 
(2004)  that  allows  the  ground  and  building  materials  to  both  reflect  and 
absorb  acoustic  energy  as  real  geologic  and  construction  materials  do. 

6.  The  porous  models  help  the  efficiency  of  the  acoustic  calculations  because 
they  produce  materials  that  do  not  control  the  time  step,  i.e.,  they  do  not 
make  the  time  step  for  convergence  smaller  and  therefore  the  analysis 
more  computationally  demanding.  Because  acoustic  propagation  in  the 
porous  materials  is  not  required  to  be  modeled  accurately,  the  number  of 
nodes  per  wavelength  based  on  the  generally  shorter  wavelengths  in  the 
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porous  materials  can  be  avoided.  The  porous  materials  use  a  flow 
resistance  that  introduces  a  non-linearity  in  the  porous-material  response. 

7.  PSTOP3D  allows  for  the  input  of  a  "geo"  model  that  defines  the 
topography  of  the  problem.  Also,  auxiliary  data  files  contain  information 
about  objects  placed  on  the  topography  such  as  buildings  or  roads. 

8.  Propagation  through  air  in  the  low  frequencies  in  the  acoustic  models  does 
not  require  a  material-attenuation  model. 

9.  The  errors  associated  with  the  second-order  solutions  for  acoustic 
propagation  are  small  compared  to  the  relative  uncertainty  of  the  ground 
or  air  properties. 

10.  Output  from  the  acoustic  analyses  is  pressures  on  2-D  slices  or  3-D 
volumes.  The  output  can  be  decimated  in  time  and  space. 

The  PSTOP3D  code  uses  data  files  to  define  the: 

1.  Grid 

2.  Topography 

3.  Propagation  medium  and  material  properties 

4.  Buildings,  roads,  pavements,  trees,  water,  etc.  (auxiliary  files) 

5.  Source  inputs 

7.1.2  Implicit/Explicit  formulation 

Explicit  and  implicit  methods  are  used  in  numerical  analysis  to  obtain 
numerical  solutions  of  time-dependent  ordinary  and  partial  differential 
equations.  PSTOP3D  uses  an  explicit  second-order-accurate,  forward- 
finite-difference  algorithm. 

Explicit  methods  use  the  state  of  the  system  at  the  current  time  to  calculate 
the  state  of  the  system  at  a  later  time.  Implicit  methods  find  a  solution  by 
solving  an  equation  involving  both  the  current  and  future  state  of  the 
system.  Implicit  methods  require  more  computations  and  can  be  harder  to 
implement,  but  they  can  use  larger  time  steps  than  an  explicit  method. 
Explicit  methods  have  the  disadvantage  of  requiring  smaller  time  steps  for 
convergence  of  the  solution.  The  determination  of  the  time  step  to  use  for 
PSTOP3D  analysis  is  critical  and  is  discussed  in  Sections  7.2  and  7.3. 

7.1.3  Sound  propagation  in  porous  medium 

From  Wilson  and  Liu  (2004),  most  common  outdoor  ground  surfaces 
cannot  be  modeled  satisfactorily  as  ideal,  rigid  surfaces.  This  is  because 
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sound  energy  propagates  into  the  pores  of  the  ground,  where  it  is 
dissipated  by  viscosity  and  thermal  conduction.  Ground  surfaces  with 
relatively  large  open  pores,  such  as  snow,  absorb  much  of  the  sound- 
energy  incident  upon  them.  Surfaces  with  small  pores,  such  as  cement  and 
asphalt,  reflect  most  of  the  energy.  The  acoustical  behavior  of  soils  is 
intermediate  between  these  extremes.  Typical  values  of  the  material 
properties  for  porous  materials  are  given  in  Table  7.1.  The  material 
properties  for  a  porous  material  needed  for  PSTOP3D  consist  of  the 
density,  wave  speed,  flow  resistivity,  porosity,  and  tortuosity. 


Table  7.1.  Typical  values  of  static-flow  resistivity,  porosity,  and  tortuosity  for 
common  porous  ground  surfaces  (Wilson  and  Liu  2004). 


Material 

Flow  Resistivity,  0 
Pa-s/m2 

Porosity,  Q 

Tortuosity,  q 

Asphalt 

3.00E+07 

0.1 

3.2 

Grass 

2.00E+05 

0.5 

1.4 

Forest 

l.OOE+05 

0.6 

1.3 

Sand 

5.00E+04 

0.35 

1.6 

Snow 

l.OOE+03 

0.6 

1.7 

The  governing  equations  for  the  acoustical  analysis  are  shown  in 
Equations  7.1  and  7.2  (Cudney  et  al.  2007,  2008). 


9p 

at 


-K,Vv  +  K,Q 


(7.1) 


dw  1 

—  =  — (Vp  +  ov)  (7.2) 

Pe 


where 

p  =  acoustic  pressure 
V  =  vector  of  particle  velocity  components 
t  =  time  (seconds) 

Ke  =  effective  bulk  modulus  of  material 
Pe  =  effective  density  of  medium  (kg/m3) 
o  =  static  flow  resistivity  (Pa/m2/s) 

=  o  for  air  (nonporous),  not  equal  to  zero  for  porous  material 
Q  =  dilation-rate  source,  radially  oscillating  sphere  (m3/sec/m3) 
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The  effective  bulk  modulus  and  effective  density  are  given  in  Equations  7.3 
and  7.4,  respectively  (Wilson  and  Liu  2004;  Cudney  et  al.  2007;  Cudney  et 
al.  2008;  Attenborough  1983). 


K. 


pc 

(7.3) 

yH 

4pq" 

3  n 

(7.4) 

where 

=  void  fraction  or  material  porosity  (volume  of  pores/unit 
volume,  m3/m3) 

q  =  tortuosity,  reflects  the  geometry  of  the  pores  (unit-less) 

Y  =  ratio  of  specific  heats  of  fluid  portion  of  porous  material, 
equals  1.4 

p  =  density  of  fluid  portion  of  porous  material  (kg/m3) 
c  =  wave  speed  in  the  fluid  portion  of  the  porous  material  (m/sec) 

The  effective  wave  speed  of  the  porous  medium  is  given  in  Equation  7.5. 


For  grass,  inserting  the  values  given  in  Table  7.1  into  Equations  7.3 
through  7.5  with  a  wave  speed  of  330  m/sec  yields: 


K. 


Pe 


pc^  _  1.2(330)"  _ 
1.4(0.5) 

4p(7"  _4l.2(l.4f 
3  n  “3  0.5 

\^e  _  /90124.1 
6.272 


90124.1Pa/m" 


6.272kg/m^ 


119.87  m  I  sec 


For  asphalt,  inserting  the  values  given  in  Table  7.1  into  Equations  7.3 
through  7.5  with  a  wave  speed  of  330  m/sec  yields: 
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K. 


pc^  1.2(330)" 


yd 

1.4(0.1) 

4pq" 

_  41.2(3.2)'  _ 

^  3  Q 

3  0.1 

1  933428.6 

ipe 

V  6.2163.8472 

163.84  kg  /  rrr 


75.48m/ sec 


As  seen  from  these  calculations,  the  effect  of  the  porous  material  is  to 
reduce  the  effective  wave  speed.  The  wave  speed  for  a  porous  material  is 
less  than  the  wave  speed  in  air.  Also,  the  wavelengths  in  the  porous 
materials  will  be  shorter  than  the  wavelengths  in  air. 

7.2  Time  Step  for  Convergence  of  Analysis 

The  stability  of  explicit  time-marching  numerical  procedures  to  solve 
partial  differential  equations  (usually  hyperbolic  partial  differential 
equations)  is  subject  to  the  Courant  time  condition.  The  numerical 
procedure  should  use  a  time  step  less  than  or  equal  to  the  Courant  time 
condition  to  insure  the  solution  converges  to  a  correct  result.  The  Courant 
time  condition  is  based  on  the  smallest  grid  increment  and  the  largest 
velocity  used  in  the  analysis.  Equation  7.6  defines  the  Courant  time  step 
(Yee  1966;  Zheng  et  al.  1999). 


At< 


1 


1 

w 


1 

A? 


(7.6) 


where 

Cmax  =  the  maximum  wave  speed  within  the  model 
Ax  =  the  minimum  grid  spacing  in  the  x  direction 
Ay  =  the  minimum  grid  spacing  in  the  y  direction 
Az  =  the  minimum  grid  spacing  in  the  z  direction 

Since  wave  speeds  are  smaller  in  porous  materials,  as  shown  in  Section 
7.1.3,  they  do  not  increase  the  time  step  required  for  stability  as  computed 
from  Equation  7.6.  However,  the  use  of  porous  materials  causes  instability 
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to  occur  as  discussed  in  Cudney  et  al.  (2007,  2008).  For  the  solution  to  be 
stable,  the  following  condition  must  be  met: 


1 


(7.7) 


This  can  be  met  by  requiring  the  time  step  to  be  less  than  or  equal  to 

But  for  a  material  such  as  asphalt  with  a  large  flow  resistivity,  the  time 
step  would  become  unreasonably  small  increasing  the  computational  cost 
by  orders  of  magnitude  for  the  analysis  of  wave  propagation  over  long 
distances.  For  asphalt,  the  time  step  would  need  to  be  below: 


At<  —  < 
a 


163.84 

3x10" 


<5x10  ^  sec 


For  a  10  sec  analysis,  asphalt  would  require  2  million  time  steps.  The 
solution  presented  by  Cudney  el  al.  (2007,  2008)  was  to  artificially  cap  the 
static  flow  resistivity  to  satisfy 


— a<l  (7.8) 

Pe 

This  was  found  to  give  good  results  for  asphalt,  which  is  the  most  reflective 
of  the  materials  listed  in  Table  7.1.  The  time  step  used  in  Equation  7.8 
would  be  smaller  than  the  Courant  time  step  but  larger  than  time  step 
computed  using  the  actual  flow  resistivity. 

The  procedure  to  cap  the  static-flow  resistivity  is  automatically  included  in 
the  PSTOP3D  program. 

7.3  Determination  of  time  step  and  SNAP  fiie  intervai 

The  time  step  used  in  PSTOP3D  is  dependent  upon  two  considerations: 

1.  The  spatial  increment  steps  must  be  small  enough  in  comparison  with  the 
wavelength  of  interest  in  order  to  make  the  numerical-dispersion  error 
negligible.  This  is  satisfied  if  there  are  10  nodes  per  wavelength. 

2.  The  time  step  must  be  small  enough  to  satisfy  the  Courant  stability 
condition  in  Equation  7.6,  which  depends  on  the  spacing  of  the  grid. 
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Assuming  lo  nodes  per  wavelength  and  a  given  grid  spacing,  the  shortest 
wavelength  that  can  be  reproduced  is: 


A„,  =10(As^,J 


(7.9) 


where  ASmin  is  the  smallest  grid  spacing  in  any  direction. 

The  maximum  frequency  that  can  reliably  be  represented  in  the  analysis  is 
dependent  upon  the  minimum  wave  speed,  Cmin,  and  the  maximum  spatial 
increment,  ASmax,  as  shown  in  Equation  7.10. 


f 

J  m 


C  C 

min  _ _ min _ 

Kax 


(7.10) 


The  minimum  wave  speed  and  maximum  grid  spacing  of  the  entire  grid 
are  used  to  compute  the  maximum  frequency  for  the  analysis.  That  is,  the 
values  produce  a  minimum  value  for  the  maximum  frequency  that  can  be 
reliably  reproduced.  Smaller  grid  spacings  or  larger  wave  speeds  will  allow 
larger  maximum  frequencies  to  be  used.  The  maximum  frequency  defined 
in  Equation  7.10  should  be  greater  than  or  equal  to  the  maximum 
frequency  of  interest  (fi).  Once  a  grid  spacing  is  determined  that  supports 
the  maximum  frequency  of  interest,  the  time  step  for  stability  can  be 
computed  using  Equation  7.6. 


A  time  increment  used  to  decimate  the  output  information  is  contained  in 
the  COMMAND. FIL  file  discussed  in  Section  7.8.  This  time  increment  is 
specified  as  a  multiple  of  the  time  step  input  for  stability  (e.g.,  output 
information  at  every  50  time  steps).  The  snap  time  interval  defines  a 
sampling  frequency  called  the  Nyquist  rate  defined  as: 


Ir 


1 

At  X  INTERVAL  SNAP 


(7.11) 


The  frequency  that  can  reliably  be  reproduced  from  this  sampling  rate  is 
called  the  Nyquist  frequency  and  is  defined  as: 


/jv  = 


Ir 

2 


(7.12) 
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The  Nyquist  frequency  for  the  snap  interval  should  be  greater  than  or 
equal  to  the  maximum  frequency  of  interest  as  shown  in  Equation  7.13. 

In  >  fi  (7-13) 

Substituting  Equation  7.11  into  7.12  and  inserting  the  result  into  Equation 
7.13  yields  Equation  7.14,  which  defines  the  snap  interval  required  to 
represent  the  maximum  frequency  of  interest. 

1 

Snap  Interval  < -  (7.14) 

^  -2xAtxfj 

7.4  Selection  of  grid  spacing  and  number  of  processors  for  PSTOP 

An  High  Performance  Computing  (HPC)  machine  contains  a  certain 
number  of  processors  per  node.  Nodes  can  be  requested  through  the  batch 
system  as  described  in  Section  7.10.  If  the  number  of  processors  requested 
is  not  a  multiple  of  the  number  of  processors  per  node,  there  will  be  idle 
processors  on  a  node.  For  example,  if  five  processors  are  requested  on  a 
machine  that  contains  four  processors  per  node,  then  two  nodes  will  be 
allocated  for  a  total  of  eight  processors.  Three  of  the  processors  will  not  be 
used.  Therefore,  to  maximize  efficiency  on  the  HPC  machine,  full  nodes 
should  be  utilized. 

PSTOP3D  requires  the  number  of  grid  points  in  the  x,  y,  and  z  directions 
to  be  evenly  divisible  by  the  number  of  processors  (also  referred  to  as  cores 
or  CPUs)  in  the  same  direction.  For  example,  for  a  grid  with  2496  nodes  in 
the  X  direction  and  16  processors  in  the  x  direction,  the  ratio  of  grid  points 
to  processors  would  be  156. 

Therefore,  the  number  of  grid  points  used  in  the  x,  y,  and  z  directions 
should  be  divisible  by  the  number  of  processors  in  the  same  direction  and 
the  total  number  of  processors  should  be  divisible  by  the  number  of 
processors  per  node  on  the  machine. 

For  example, 

1.  Assume  that  the  grid  is  2336  x  2328  x  508  nodes  in  x,  y,  z  directions. 

2.  Assume  32  processors  per  node. 

3.  Use  16  X 12  X  4  processors  in  the  x,  y,  z  directions. 


ERDC  TR-15-5 


114 


4.  X  direction,  2336  / 16  =  146  FD  nodes/processor 

5.  y  direction,  2328/12  =  194  FD  nodes/processor 

6.  z  direction,  508/4  =  127  FD  nodes/processor 

7.  Total  number  of  processors  =  16  x  12  x  4  =  786  processors 

8.  (Total  number  of  processors)/ (processors  per  node)  =  786/32  =  24 

9.  Request  a  total  of  768  processors  (24  nodes). 

10.  In  the  pstop3d.h  file  discussed  in  Section  7.6,  the  number  of  processors  in 
the  X,  y,  and  z  directions  would  be  16  x  12  x  4. 

PSTOP3D  also  has  a  restriction  that  must  be  observed  concerning  the 
number  of  processors  in  the  z  direction.  PSTOP3D  requires  that  the  PML 
region  must  be  contained  within  a  single  processor.  Also,  for  non-flat 
topography,  the  parallel-to-topography  output  must  be  contained  within  a 
single  processor.  PSTOP3D  will  alert  the  user  if  these  conditions  are  not 
met  and  halt  execution.  If  this  occurs,  the  number  of  processors  used  in 
the  z  direction  would  have  to  be  reduced. 

7.5  Grid  spacing,  time  step,  and  snap  intervai  for  urban  area 

The  modeling  of  the  urban  area  described  in  this  report  focused  on  the 
propagation  effects  of  infrasound  signals  observed  in  the  urban  environ¬ 
ment.  Infrasound  is  sound  that  is  20  Hz  and  below.  For  the  analyses 
performed,  the  maximum  frequency  of  interest  was  chosen  to  be  15  Hz. 

The  grid  size,  time  step,  and  snap  interval  were  determined  by  the 
following  procedure: 

1.  The  modeling  of  buildings  was  an  important  aspect  of  the  analysis,  so  a 
spacing  of  approximately  1  m  was  selected  as  a  maximum  spacing  to  use  to 
adequately  reflect  the  geometry  of  the  buildings. 

2.  fmax  was  computed  from  Equation  7.10  using  a  minimum  wave  speed  of 
330  m/sec  and  a  grid  spacing  of  1  m.  This  frequency  (32.9  Hz)  was  greater 
than  the  maximum  frequency  of  interest  (15  Hz). 

3.  The  time  step  required  for  stability  (i.e.,  the  Courant  time  step)  was 
computed  using  Equation  7.6.  The  maximum  wave  speed  and  minimum 
grid  spacing  in  the  x,  y,  and  z  directions  were  used  to  compute  the  time 
step  of  0.00175  sec. 

4.  Due  to  the  use  of  porous  materials  in  the  model,  the  time  step  was 
required  to  be  less  than  5  x  10-^  sec  as  discussed  in  Section  7.2.  The  use  of 
this  time  step  was  deemed  to  be  unfeasible.  Therefore,  a  larger  time  step 
was  selected  and  the  stability  of  the  algorithms  was  assured  by  letting 
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PSTOP3D  cap  the  flow  resistivities  as  discussed  in  Section  7.2.  A  time  step 
of  0.0005  sec  was  selected  as  a  reasonable  time  step  and  is  less  than  the 
Courant  time  condition.  This  resulted  in  20,500  time  steps  being  run  to 
allow  the  waves  to  propagate  across  the  complete  urban  area. 

5.  The  snap  interval  was  computed  using  Equation  7.14  using  the  time  step 
of  0.0005  sec  and  the  maximum  frequency  of  interest  of  15  Hz.  This 
resulted  in  a  maximum  snap  interval  of  66.7.  Therefore,  a  snap  interval  of 
50  was  used  in  the  analyses. 

6.  Based  on  the  grid,  a  frequency  of  32.9  Hz  could  be  represented.  Based  on 
the  analysis  time  step,  a  frequency  of  2000  Hz  could  be  represented.  Based 
on  the  snap  interval,  a  frequency  of  20  Hz  could  be  represented.  Therefore, 
the  maximum  frequency  for  the  analysis  that  could  be  represented  was 

20  Hz. 

7.6  Edit  the  pstopSd.h  file  for  PST0P3D 

To  compile  PSTOP3D,  a  header  file  is  required  for  the  FORTRAN 
program.  This  header  file  contains  information  about  the  dimensions  of 
arrays  specific  to  the  problem  being  analyzed.  This  header  file  is  created  by 
the  MATLAB  script  rvloadset_io.m  as  discussed  in  Section  6.6.  This 
header  file  must  be  edited  for  the  number  of  processors  used  in  the 
analysis  as  discussed  in  Section  7.4. 

The  header  file  has  a  name  like  PSTOPsD.h. RUNID.txt  where  RUNID 
has  the  form  MODELID_LOADID.  MODELID  is  based  on  the  model 
information  contained  in  the  spreadsheet  discussed  in  Section  5.3  and 
LOADID  is  taken  from  the  analysis  identifier  entered  when  running  the 
script  rvloadsdet_io.m  as  discussed  in  Section  6.6.  The  MODELID  is  a 
concatenation  of  the  id.model,  id.grd,  id.topo,  and  id.med  data  items. 
The  file  name  shown  in  Section  6.6  was  pstop3d.h.mctp_URi5.txt 
where  the  MODELID  is  mctp  and  the  LOADID  is  UR15.  The  header  file 
is  shown  below  in  Courier  font.  Green  text  denotes  text  that  must  be 
edited  and  red  text  is  additional  explanation. 


C  ===  START  OF  FILE:  pstopSd.h,  written  by  Matlab  m-file  rvloadset_V.m; 
C  ===  edit  where  indicated  to  correct  undesired  default  values. 

C  S.A.  KETCHAM 

C  USA  ERDC  CRREL,  72  LYME  RD,  HANOVER  NH  03755 
C  EMAIL : STEPHEN . A . KETCHAM@USACE . ARMY . MIL 
C  VOICE:  603-646-4601 

C  PURPOSE: 

C  INCLUDE  FILE,  DEFINES  PARAMETERS  FOR  PSTOP3D_*.F 
C  PARALLEL  MPI  FD-TD  VISCOELASTIC  OR  ACOUSTIC  PROPAGATION  MODEL 
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C  ****GENERAL  DEFINITIONS**** 

C  nEW_mesh  =  NUMBER  OF  PROCESSORS  IN  X  (EAST-WEST)  DIRECTION 
C  nNS_mesh  =  NUMBER  OF  PROCESSORS  IN  Y  (NORTH-SOUTH)  DIRECTION 
C  nUD_mesh  =  NUMBER  OF  PROCESSORS  IN  Z  (UP-DOWN)  DIRECTION 

C  NNX  =  NUMBER  OF  SINGLE-PROCESSOR-DOMAIN  GRID  POINTS  IN  X  (EAST-WEST)  DIRECTION 
C  NNY  =  NUMBER  OF  SINGLE-PROCESSOR-DOMAIN  GRID  POINTS  IN  Y  (NORTH-SOUTH)  DIRECTION 
C  NNZ  =  NUMBER  OF  SINGLE-PROCESSOR-DOMAIN  GRID  POINTS  IN  Z  (UP-DOWN)  DIRECTION 
C  LO  =  LENGTH  OF  THE  OVERLAP  BETWEEN  PROCESSOR  DOMAINS  FOR  SPATIAL  DIFFERENCES 
C  NX_MODEL_SIZE+2*LO  =  NUMBER  OF  FULL-MODEL  GRID  POINTS  IN  X  (EAST-WEST)  DIRECTION 
C  NY_MODEL_SIZE+2*LO  =  NUMBER  OF  FULL-MODEL  GRID  POINTS  IN  Y  (NORTH-SOUTH) 
DIRECTION 

C  NZ_MODEL_SIZE+2*LO  =  NUMBER  OF  FULL-MODEL  GRID  POINTS  IN  Z  (UP-DOWN)  DIRECTION 
C  NREL  =  NUMBER  OF  VISCOELASTIC  RELAXATION  MECHANISMS 

C  NNP  =  WIDTHS  OF  PML  ABSORBING  ANULUS  ALL  MODEL  SIDES  (<=NNX, NNY, NNZ ) ; 

C  NRANK  =  RANK  OF  STRESS  (OR  PRESSURE)  TENSOR  TO  DIMENSION  PML  MATRIX  SIZE; 

C  ================================================================= 

Q  -k-k-k-k-k  important  !  !  !  !  ***** 

C  NX_MODEL_SIZE  /  NNX  MUST  EQUAL  A  WHOLE  NUMBER  BEFORE  ROUNDING! 

C  NY_MODEL_SIZE  /  NNY  MUST  EQUAL  A  WHOLE  NUMBER  BEFORE  ROUNDING! 

C  NZ_MODEL_SIZE  /  NNZ  MUST  EQUAL  A  WHOLE  NUMBER  BEFORE  ROUNDING! 

C  E.G. ,  THIS  WORKS: 

C  NX_MODEL_SIZE=I20 
C  nEW_MESH=3 

C  NX_MODEL_SIZE/nEW_MESH=40=NNX 
C  THIS  DOES  NOT  WORK: 

C  NX_MODEL_SIZE=121 
C  nEW_MESH=3 

C  NX_MODEL_SIZE/nEW_MESH=40=NNX 

C  ================================================================= 


C  EDIT  THE  PARAMETER  STATEMENTS  BELOW  TO  DESIRED  ANALYSIS  SPECIFICATIONS 


C  SPECIFY  PROCESSOR  MESH 
PARAMETER  (nEW_mesh=R) 
PARAMETER  (nNS_mesh=I d) 
PARAMETER  (nUD  mesh=^  ! 


(DEFAULT  IS  4 -CORE  DOMAIN) : 
!  SIZE  IN  X-DIRECTION,  i.e. 
!  SIZE  IN  Y-DIRECTION,  i.e. 
SIZE  IN  Z-DIRECTION,  i.e., 


EAST-WEST  DIMENSION 
NORTH-SOUTH  DIMENSION 
UP-DOWN  DIMENSION 


The  values  of  the  size  of  the  processors  in  the  x,  y,  and  z  directions  should  be 
determined  as  discussed  in  Section  7.4.  This  file  was  composed  for  a  32 
processor/node  machine. 


C  SPECIFY  WIDTH  OF  PML  ANULUS  SURROUNDING  MODEL  ALL  SIDES: 
PARAMETER  (NNP=10) 


This  value  of  NNP  should  remain  10.  A  value  of  10  for  the  PML  region  is  hardcoded 
in  the  MATLAB  routines.  This  is  discussed  in  Section  5.6.  The  PML  region  must  be 
contained  within  a  single  processor  level. 

C  SPECIFY  NUMBER  OF  VISCOELASTIC  MECHANISMS  (CHANGE  DEFAULT  IN  ONE  CIRCUMSTANCE 
C  ONLY:  TO  PERFORM  ELASTIC  ANALYSIS  WITH  VISCOELASTIC  MODEL,  SET  NREL=0) : 
PARAMETER  (NREL=0) 


C  ================================================================= 

C  DO  NOT  EDIT  THE  PARAMETER  STATEMENTS  BELOW 
C  ENTIRE  MODEL'S  GRID  SIZE: 

PARAMETER  (NX_MODEL_SI ZE=2 4 9 6 )  This  is  the  number  of  x  nodes  in  the  mesh. 
PARAMETER  (NY_MODEL_SI ZE=2 4 9 6 )  This  is  the  number  of  y  nodes  in  the  mesh. 
PARAMETER  (NZ_MODEL_SIZE=256)  This  is  the  number  of  z  nodes  in  the  mesh. 

If  the  topography  is  not  flat,  the  parallel-to-topography  output  must  be  contained 
within  a  single  z-processor  level. 

C  I_SEISMIC_FLAG=1 :  SEISMIC  ANALYSIS;  I_SEISMIC_FLAG=0 :  ACOUSTIC  ANALYSIS 
PARAMETER  (I  SEISMIC  FLAG=0) 
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C  DOMAIN  GRID  SIZE  FOR  EACH  PROCESSOR  IS  CALCULATED — DO  NOT  CHANGE: 
PARAMETER  (NNX=NX_MODEL_SIZE/nEW_mesh) 

PARAMETER  (NNY=NY_MODEL_SIZE/nNS_mesh) 

PARAMETER  (NNZ=NZ_MODEL_SIZE/nUD_mesh) 

C  OVERLAPS: 

PARAMETER  (LO=l) 

PARAMETER  (LOA=I) 

PARAMETER  (LOE=0) 

PARAMETER  (LOV=0) 

C  VARIABLE-DIMENSION  PARAMETERS: 

P ARAME  TER  ( NNX A=NNX , NN Y A=NN Y , NN  Z A=NN  Z ) 

PARAMETER  (NNXE=1 , NNYE=1 , NNZE=1 ) 

PARAMETER  (NNXV=1 , NNYV=1 , NNZV=1 ) 

PARAMETER  (NREL_DIMPARAM=1 ) 

PARAMETER  (NNPA=NNP, NNPE=1 , NNPV=1 ) 

C  RANK  OF  STRESS  TENSOR  (=2  FOR  STRESS,  =0  FOR  PRESSURE) : 

PARAMETER  (NRANK=0) 

C  MATERIAL  INDEX  DEFINITION  TYPE  CAN  BE,  E.G.,  BY  NODES,  LAYERS,  ETC., 
C  AS  DEFINITION  CAPABILITIES  EVOLVE--DO  NOT  CHANGE: 

PARAMETER  (MATL_DEF_TYPE=1 ) 

C  NUMBER  OF  LAYERS  TO  DEFINE  MATERIAL  INDICES — DO  NOT  CHANGE: 

PARAMETER  (MATL_DEF_PARAMETER_1=4 ) 


C  ================================================================= 

C  EDIT  THE  PARAMETER  STATEMENTS  BELOW:  SET  TO  1  TO  SAVE  MEMORY  WHEN  NOT 
CALCULATING  STRAIN 

C  VARIABLE-DIMENSION  PARAMETERS  FOR  STRAIN  QUANTITY  CALCULATION 
PARAMETER  (NNXEE=1 , NNYEE=1 , NNZEE=1 ) 


C  ================================================================= 

C  DO  NOT  EDIT  PARAMETER  STATEMENTS  BELOW 

C  LOAD  NUMBER  I:  AI .  NOMINAL  PRESSURE,  rvloadset_window_FRicker 
C  NUMBER  OF  APPLIED  SOURCES  (I.E.,  BODY  FORCES,  VELOCITY  CONSTRAINTS,  OR  Q 
SOURCES) =I 

C  NUMBER  OF  g(t)  TIME  STEPS=299 

C  SOURCE  POSITION  (xs , ys , height ) = ( 30 . 0 4 8 I , 2 4 6 9 . 95 , 1 . 4 034 9 ) ;  height  IS  HEIGHT  ABOVE 
TOPOGRAPHIC  SURFACE 

C  SOURCE  POSITION  (xs, ys, zs) = (30 . 0481, 2469 . 95, 38) ;  zs  IS  Z-COORDINATE  POSITION  AND 
ELEVATION  ABOVE  DATUM  =  z  =  0  m 

C  VECTOR  INPUT  FLAG  (FLAG=I  FOR  VECTOR  INPUT,  0  FOR  SCALAR  INPUT) 

PARAMETER  ( I_VECTOR_INP_FLAG=0 ) 

C  NUMBER  OF  LOADS  IN  LOAD  SET 
PARAMETER  (NUM_LOADS=I ) 

C  TOTAL  NUMBER  OF  SOURCES  IN  ALL  LOADS  OF  LOAD  SET 
PARAMETER  (NUM_FORCES=I ) 

C  TOTAL  NUMBER  OF  G(t)  STEPS  IN  ALL  LOADS  OF  LOAD  SET 
PARAMETER  (NUM_GSTEPS=2 9 9 ) 

The  NUM_GSTEPS  is  for  the  specific  source.  If  the  source  changes,  this  number  will 
change.  The  files  defining  the  source  are  written  by  the  rvloadset.m  script 
described  in  Section  6.6.  The  COMMAND. FIL  file  described  in  Section  7.8  lists  the 
files  for  use  in  PSTOP3D. 


7.7  Compile  PSTOP3D 

To  compile  PSTOP3D,  the  pstop.h.RUNID.txt  file  should  be  copied  to 
the  directory  where  the  PSTOP  Fortran  files  are  located  on  the  HPC 
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machine.  The  RUNID  part  of  the  file  name  is  explained  in  Section  7.6. 
The  following  script  can  be  run  to  compile  PSTOP3D: 


# ! /bin/ksh 

#  pstop3d  compile  for  Garnet  script 

#  754=rwxr-xr--  Group  can  read  and  execute 

#  Use  "chmod  754  pstop3d_compile . sh"  if  not  executable 

#  Current  directory  must  have  this  script,  pstop3d . h . $1 . txt ,  and  all  pstop3d  .f 
files 

#  Current  directory  will  be  the  output  directory 

#  $1  is  the  analysis  and  pstop3d  executable  ID  (long  form,  i.e.,  not  necessarily 
the  . SNP  file  prefix) 

#  To  submit  for  analysis  mctp_UR15  (example  below) ,  with  input  for  $1 

#  pstop3d_compile . sh  mctp_UR15 

#  Initiate  log  entries 

date  2>&1  |  tee  -a  pstop3d_compile_$l . log 
echo  "using  file  pstop3d . h . $1 . txt" 


#  Copy  include  .txt  file  for  current  analysis 

cp  pstop3d. h . $1 . txt  pstop3d.h  2>&1  |  tee  -a  pstop3d_compile_$l . log 


#  Compile  code 
echo  "compiling. . . " 
ftn  -fast  -byteswapio 
ftn  -fast  -byteswapio 
ftn  -fast  -byteswapio 
ftn  -fast  -byteswapio 
ftn  -fast  -byteswapio 
ftn  -fast  -byteswapio 


-c  pstop3d_main . f  2>&1  |  tee  -a  pstop3d_compile_$l . log 
-c  pstop3d_mpi . f  2>&1  |  tee  -a  pstop3d_compile_$l . log 

-c  pstop3d_abc . f  2>&1  |  tee  -a  pstop3d_compile_$l . log 

-c  pstop3d_strs . f  2>&1  |  tee  -a  pstop3d_compile_$l . log 
-c  pstop3d_pak . f  2>&1  |  tee  -a  pstop3d_compile_$l . log 

-c  pstop3d_vel . f  2>&1  |  tee  -a  pstop3d_compile_$l . log 


echo  "linking. . . " 

ftn  -fast  -byteswapio  -o  pstop3d_$l  *.o  2>&1  |  tee  -a  pstop3d_compile_$l . log 


The  script  may  be  executed  by  typing: 

./compile.ksh  RUNID 

where  RUNID  is  the  designator  for  the  analysis  defined  by  the 
rvloadset_io.m  script  (e.g.,  mctp_URi5). 

This  script  does  the  following: 

1.  Copies  pstop3d.h.RUNID.txt  to  pstop.h. 

2.  Compiles  the  FORTRAN  modules  and  links  them  creating  the  executable. 

3.  The  executable  file  pstop3d_RUNID  is  made,  e.g., 

pstop3d_mctp_URi5 . 

4.  This  executable  is  “good”  for  the  grid  size  and  load  set  defined.  Therefore, 
if  all  load  sets  are  the  same  and  just  vaiy  in  their  location,  then  this  file  is 
the  same.  The  executable  will  be  the  same  and  can  be  used  for  multiple 
loadings. 
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7.8  Edit  the  COMMAND.FIL  file  for  PST0P3D 

One  of  the  files  created  by  the  rvloadset_io.m  script  is  a  file  named 
COMMAND.FIL.RUNID.txt  where  the  RUNID  is  the  run  identifier 
described  in  Section  7.6.  This  file  controls  the  run  parameters  for 
PSTOP3D  and  needs  to  be  edited  for  the  specific  problem  being  analyzed. 
The  file  is  written  for  a  specific  load  set  as  defined  by  the  rvloadset_io.m 
script.  The  file  is  shown  below  in  Courier  font.  Green  text  denotes  text 
that  must  be  edited  while  red  text  is  additional  explanation. 


Parameter  inputfile  COMMAND.FIL  for  pstopSd  written  by  Matlab  m-file 
rvloadset  V.m. 


Ittctp  UR15|  ;  SNAP_FILE.  EDIT  TO  DISTINGUISH  OUTPUT  FILES  IF  NEEDED. 

<15  CHARACTERS. 


This  will  be  the  name  of  the  SNAP  files  created  during  the  analysis.  The  SNAP 
files  are  files  that  contain  snapshots  of  the  output  at  selected  time  increments. 


mctp_geo . dat 
1.0016,1.0016,1 
ONLY. 

0.0005 


GEO_EILE.  DO  NOT  EDIT. 

Dxi, Dkappa, Deta  (m) .  DO  NOT  EDIT:  FOR  INFO-FILE  OUTPUT 
DT  (s) .  DO  NOT  EDIT. 


This  is  the  time  step  for  the  analysis  which  must  satisfy  the  Courant  condition 
discussed  in  Section  7.2. 

|20500|  ;  NT.  EDIT  TO  SET/CHANGE  DURATION  OE  ANALYSIS. 

This  is  the  number  of  time  steps  in  the  analysis. 


0 


OF  OUTPUT. 


;  Blank  spot  for  future  pml  parameter. 

;  INTERVAL_SNAP .  EDIT  TO  ENSURE  DESIRED  NYQUIST  EREQUENCY 


The  value  of  the  SNAP  interval  should  be  selected  to  insure  the  desired  Nyquist 
frequency  as  described  in  Section  7.3.  This  is  a  multiple  of  DT .  Therefore,  the 
time  increment  of  the  output  for  this  case  is  DT*INTERVAL_SNAP . 

i.626, 545, 31|  ;  IX_SNAP_PLN,  IY_SNAP_PLN,  IZ_SNAP_PLN .  EDIT  TO  SELECT  CD- 

OUTPUT  PLANES  (0=NO  SAVE) . 

The  IX_SNAP_PLN,  IY_SNAP_PLN,  IZ_SNAP_PLN  values  are  the  nodes  through  the  point 
of  application  of  the  source.  This  can  be  changed  to  select  another  location. 

1,1,1  ;  NDX_SNAP_2D,NDY_SNAP_2D,NDZ_SNAP_2D:  EDIT  TO  SPATIALLY 

DECIMATE  2D-OUTPUT  PLANES. 

I  ;  I_TOPO_OUTPUT_FLAG.  SET=0  FOR  Z=Z ( IZ_SNAP_PLN)  SLICE; 

SET=I  FOR  PARALLEL-TO-TOPO  "SLICE." 

I  ;  NODES_ABOVE_TOPO.  EDIT  TO  SELECT  PARALLEL-TO-TOPO  2D 

OUTPUT.  (IGNORED  WHEN  I_TOPO_OUTPUT_ELAG  =  0) . 

The  TOPO_OUTPUT_FLAG  is  set  for  parallel  to  topography;  therefore,  this  line 
selects  the  results  at  1  node  above  the  topography.  If  the  aux . options . modifyTopo 
flag  as  discussed  in  Section  5.7  for  buildings  is  set,  the  results  are  output  1 
node  above  the  tops  of  the  buildings. 

|1, 248  61  ;  IX_BEG_SNAP_2D,  IX_END_SNAP_2D.  EDIT  TO  CHOOSE  INDICES  FOR 

BOUNDING  2D  OUTPUT. 

tn ,  2 4  8  6|  ;  IY_BEG_SNAP_2D,  IY_END_SNAP_2D.  EDIT  TO  CHOOSE  INDICES  EOR 

BOUNDING  2D  OUTPUT. 

|1, 2  4  6|  ;  IZ_BEG_SNAP_2D,  IZ_END_SNAP_2D.  EDIT  TO  CHOOSE  INDICES  FOR 

BOUNDING  2D  OUTPUT. 
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The  node  numbers  for  the  2D  output  should  be  adjusted  to  accommodate  the  10  node 
PML.  So  for  2496  nodes,  begin  at  11  and  stop  at  2486.  For  256  nodes,  begin  at  11 
and  end  at  246. 


1,2496 

BOUNDING  3D  OUTPUT. 
1,2496 

BOUNDING  3D  OUTPUT. 
1,256 

BOUNDING  3D  OUTPUT. 


IX_BEG_SNAP_3D, IX_END_SNAP_3D.  EDIT  TO  CHOOSE  INDICES  EOR 
IY_BEG_SNAP_3D, IY_END_SNAP_3D.  EDIT  TO  CHOOSE  INDICES  EOR 
IZ_BEG_SNAP_3D, IZ_END_SNAP_3D.  EDIT  TO  CHOOSE  INDICES  EOR 


0,0,0 

DECIMATE  3D  OUTPUT 


;  NDX_SNAP_3D,NDY_SNAP_3D,NDZ_SNAP_3D:  EDIT  TO  SPATIALLY 
(any  0=NO  SAVE) . 


0,0, 0,1  ;  Indicators  for  saving  quantities  U,V,W,P.  (1=SAVE,  0=NO 

SAVE) . 

|ttctp~UR15  Is.daj  ;  SOURCESET_FILE .  DO  NOT  EDIT  EXCEPT  IN  SPECIAL 

CIRCUMSTANCES.  <28  CHARACTERS. 

^ctp  UR15  gt.daj  ;  GT_EILE.  DO  NOT  EDIT  EXCEPT  IN  SPECIAL  CIRCUMSTANCES.  <28 

CHARACTERS . 


The  above  file  names  are  for  the  source  defined  when  the  rvloadset_10 . m  script  is 
run.  This  can  be  changed,  but  this  also  affects  lines  in  the  pstopSd.h  header  file 
as  discussed  in  Section  7.6. 


GLOSSARY  FOR  INPUT  TERMS  ABOVE 


SNAP_FILE:  ROOT  OF  OUTPUT  FILE  NAME. 

GEO_FILE:  TOPOGRAPHY  AND  VELOCITY  MODEL  INPUT  FILE  NAME. 

Dxi, Dkappa, Deta  (m) :  MIN  X,Y,Z  SPACING  OF  THIS  ANALYSIS'  RECTANGULAR  F-D  GRID. 

DT  (s) :  TIME  STEP.  MUST  MEET  STABILITY  CONDITION.  MUST  MATCH  DT  IN  GT_FILE. 

NT:  TOTAL  NUMBER  OF  TIME  STEPS  (IF  ANALYSIS  COMPLETES) . 

Euture  pml  coefficient  if  needed. 

INTERVAL_SNAP :  TIME  DECIMATION  ORDER  FOR  SNAPSHOT  OUTPUTS. 

IX_SNAP_PLN, IY_SNAP_PLN, IZ_SNAP_PLN:  PLANE  INDICES  FOR  2D  SLICE  OUTPUT 
(0=NO  SAVE).  OUTPUT  WILL  BE  AT  X ( IX_SNAP_PLN) ,  Y ( IY_SNAP_PLN) ,  AND 
Z (IZ_SNAP_PLN) ,  WHERE  THE  COORDINATES  HERE  USE  PSTOP3D  INDICES 
THAT  START  AT  1-LO,  WHICH  ARE  NOT  TO  BE  CONFUSED  WITH  rvgeo_V.m 
INDICES  THAT  START  AT  1.  THE  CONVERSION  BETWEEN  THESE  TWO  INDEXING 
SYSTEMS  IS  index_pstop=index_rvgeo-LO.  ALSO  Z=CONSTANT  OUTPUT  WILL  BE 
OVERRIDDEN  BY  PARALLEL-TO-TOPO-SURFACE  OUTPUT  IF  I_TOPO_OUTPUT_FLAG=l . 

NDX_SNAP_2D,NDY_SNAP_2D,NDZ_SNAP_2D:  SPATIAL  DEC  ORDERS  FOR  2D  SNAPSHOTS. 

I_TOPO_OUTPUT_FLAG:  =  0  TO  OUTPUT  Z=Z ( IZ_SNAP_PLN)  SLICE;  =  1  TO  OUTPUT 
PARALLEL-TO-TOPO-SURFACE  "SLICE" . 

NODES_ABOVE_TOPO:  =  NUMBER  OF  Z-DIRECTION  NODES  BETWEEN  THE  TOPO  LAYER  AND  THE 
X-Y  OUTPUT  "SLICE"  (IGNORED  WHEN  I_TOPO_OUTPUT_FLAG  =  0) .  USE  POSITIVE 
NUMBER  TO  DESIGNATE  OUTPUT  FROM  ABOVE  TOPO  SURFACE,  AND  A  NEGATIVE  NUMBER 
TO  DESIGNATE  OUTPUT  FROM  ABOVE  TOPO  SURFACE. 

IX_BEG_SNAP_2D, IX_END_SNAP_2D:  INDICES  FOR  BOUNDING  2D  X-DIRECTION  OUTPUT 
BY  X (IX_BEG_SNAP_2D)  AND  X ( IX_END_SNAP_2D) ,  WHERE  THE  COORDINATES  HERE  USE 
PSTOP3D  INDICES  THAT  START  AT  I-LO.  THE  CONVERSION  BETWEEN  THE  PSTOP  AND 
rvgeo_V.m  INDEXING  SYSTEMS  IS  index_pstop=index_rvgeo-LO . 

IY_BEG_SNAP_2D, IY_END_SNAP_2D:  INDICES  FOR  BOUNDING  2D  Y-DIRECTION  OUTPUT. 

IZ_BEG_SNAP_2D, IY_END_SNAP_2D:  INDICES  FOR  BOUNDING  2D  Z-DIRECTION  OUTPUT. 

IX_BEG_SNAP_3D, IX_END_SNAP_3D:  INDICES  FOR  BOUNDING  3D  X-DIRECTION  OUTPUT. 

IY_BEG_SNAP_3D, IY_END_SNAP_3D:  INDICES  FOR  BOUNDING  3D  Y-DIRECTION  OUTPUT. 

IZ_BEG_SNAP_3D, IY_END_SNAP_3D:  INDICES  FOR  BOUNDING  3D  Z-DIRECTION  OUTPUT. 

NDX_SNAP_3D,NDY_SNAP_3D,NDZ_SNAP_3D:  SPATIAL  DECIMATION  ORDERS  FOR  3D  X,Y,Z 
OUTPUTS  RESPECTIVELY. 

IU_SNAP, IV_SNAP, IW_SNAP, IP_SNAP:  U,V,W,P  OUTPUT  INDICATORS  (1=SAVE,  0=NO  SAVE). 

LOADSET_FILE :  FILE  WITH  SETS  OF  BODY  FORCES,  VELOCITY  CONSTRAINTS,  OR  Q  SOURCES 

AND  LOCATIONS  FOR  EACH  LOAD. 

GT  FILE:  FILE  WITH  TIME  SERIES  GIVING  FINAL  AMPLITUDE  VARIATION  FOR  EACH  LOAD. 
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7.9  Edit  the  COMMAND.AUX  file  for  PST0P3D 

PSTOP3D  uses  files  termed  auxiliary  files.  These  files  define  auxiliary 
structures  that  are  placed  on  or  under  the  topography  such  as  buildings, 
bridges,  and  paved  areas.  The  modeling  of  auxiliary  structures  is  discussed 
in  Sections  3.3,  3.4,  and  3.5.  The  definitions  of  auxiliary  structures  in  the 
Excel  spreadsheets  are  discussed  in  Sections  5.7  and  5.8.  The  processing  of 
auxiliary  structures  using  MATLAB  is  discussed  in  Section  6.5. 

When  the  auxiliary  structures  are  processed  by  MATLAB,  a  file  named 
COMMAND  .AUX.ID_AUX.legacy.txt  is  created.  This  file  is  created  in 
a  subdirectory  named  legacy.  The  ID_AUX  part  of  the  file  name  is  a 
concatenation  of  the  MODELID  described  in  Section  5.3,  the  text  _aux_, 
and  the  id.aux  data  item.  The  ID  is  a  concatenation  of  the  id.model, 
id.grd,  id.topo,  and  id.med  data  items  (e.g.,  mctp_aux_i). 

This  file  lists  the  names  of  the  auxiliary  files  created  as  defined  in  the  Excel 
spreadsheets.  If  one  auxiliary  spreadsheet  is  processed,  the  file  will  contain 
one  file  name,  if  five  auxiliary  spreadsheets  are  processed  in  the  same 
execution  of  the  MATLAB  rvgModel.m  script  as  described  in  Section  5.3, 
then  five  file  names  will  be  written  to  the 
COMMAND.AUX.ID_AUX.legacy.txt  file. 

If  auxiliary  spreadsheets  are  processed  one  at  a  time,  the 
COMMAND .AUX.ID_AUX.legacy.txt  file  will  only  contain  the  file 
name  generated  by  the  last  auxiliary  spreadsheet  processed.  Therefore,  the 
file  will  need  to  be  edited  to  include  the  names  of  all  auxiliary  files 
generated  and  to  be  used  in  the  analysis.  Only  the  desired  auxiliary  files 
for  a  particular  analysis  should  be  included  in  this  file.  For  example,  a  set 
of  auxiliary  files  could  be  used  to  define  buildings  as  individual  buildings 
and  another  set  could  be  used  to  define  buildings  aggregated  into  larger 
areas  representing  the  buildings.  This  would  require  two  different 
COMMAND .AUX.ID_AUX.legacy.txt  files  to  describe  these  sets  of 
buildings  for  analysis. 

The  auxiliary  files  created  have  names  like  mctp_aux_1building2.dat. 
The  first  part  of  the  name  is  the  ID_AUX  identifier  discussed  previously. 
The  last  part  of  the  file  name  is  based  on  the  name  of  the  spreadsheet 
defining  the  auxiliary  structures  as  discussed  in  Section  5.3. 
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The  COMMAND .AUX.ID_AUX.legacy.txt  is  shown  below  in 
Courier  font.  Green  text  denotes  text  that  must  be  edited. 

COMMAND. AUX  provides  auxiliary  model  data  for  PSTOP3D:  #  files,  file  names, 
max  DT .  Later  data  overwrites  earlier  data--order  aux  files  accordingly. 


191 _  ;  Number  of  matl  &  indices  auxiliary-file  pairs 

Bctp_aux_l PavedA . da^  ;  Prefix  for  matl  &  indices  auxiliary  files 

Bctp_aux_lPavedB .  dal| 

Bctp_aux_lPavedC .  datj 
Bctp_aux_lPavedD .  da1| 

Bctp_aux_lPavedE .  da1| _ 

Bctp_aux_l PavedA- 1 s lands  .  dal| 

Bctp_aux_lPavedB-Islands .  dalj 
Bctp_aux_lPavedC-Islands .  dalj 
Bctp_aux_lPavedD-Islands .  dalj 
Bctp  aux  IPavedE-Islands  .  daij 
Bctp_aux_lbuildingl 
Bctp_aux_lbuilding2 
Bctp_aux_lbuilding3 
Bctp_aux_lbuilding4 
Bctp_aux_lbuilding5 
Bctp_aux_lbuilding6 
Bctp_aux_lbuilding7 
Bctp_aux_lbuilding8 
Bctp  aux  IBridge.da 


0.00764907 


=  max  time  step  (s)  for  these  files 


Modify  the  above  file  to  list  all  auxiliary  files  in  the  correct  order.  Later 
files  overwrite  previous  file  definitions.  Therefore,  the  paved  areas  are 
listed  first.  These  areas  define  larger  areas  than  needed  because  they  do 
not  account  for  islands  as  discussed  in  Section  5.8.  The  islands  are  defined 
next  to  set  portions  of  the  larger  paved  areas  defined  as  asphalt  back  to 
grass.  Auxiliary  files  defining  buildings  are  listed  next.  Finally,  the 
auxiliary  file  defining  bridges  is  listed. 

The  time  step  shown  in  the  file  corresponds  to  the  Courant  condition 
computed  for  the  auxiliary  materials.  This  value  is  for  the  material  that  the 
file  was  written.  That  is,  if  an  auxiliary  structure  is  processed  one 
spreadsheet  at  a  time,  this  is  for  the  last  spreadsheet  processed.  If  several 
spreadsheets  are  processed,  this  is  the  smallest  time  computed  for  all  of 
the  materials  specified  in  the  spreadsheets. 

7.10  How  to  execute  PST0P3D 

After  the  PSTOP3D  program  has  been  compiled  as  discussed  in  Section 
7.7,  the  data  files  that  PSTOP3D  requires  must  be  uploaded  to  the  HPC 
machine.  The  files  needed  are: 
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5.  The  load  set  files  produced  by  the  rvloadset_io.m  script  discussed  in 
Section  6.6.  These  files  are  listed  in  the  COMMAND.FIL  file  discussed  in 
Section  7.8.  These  files  are  named  RUNID_ls.dat  and  RUNID_gt.dat 
where  the  RUNID  is  defined  in  Section  7.6. 

6.  The  COMMAND.FIL.RUNID.txt  file  (e.g., 
COMMAND.FIL.mctp_UR15.txt) 

7.  The  COMMAND.AUX.ID_AUX.legacy.txt  file  (e  g., 
COMMAND.AUX.mctp_aiix_i.legacy.txt) 

8.  All  files  beginning  with  RUNID  and  ID_AUX  in  the  legacy  directory. 
RUNID  is  described  in  Section  7.6  and  ID_AUX  in  Section  7.9.  For 
example,  if  RUNID  is  mctp  and  ID_AUX  is  mctp_aiix_i  then  all  files 
starting  with  mctp_geo.dat  and  mctp_aiix_i  should  be  copied.  This  is 
all  files  containing  information  about  the  topography  and  auxiliaiy 
structures. 

The  files  should  be  copied  to  an  analysis  directory  on  the  HPC  machine. 
The  PSTOP3D  executable  should  be  copied  to  the  same  directory. 

Once  all  required  files  are  copied  to  a  directory,  a  batch  script  can  be 
submitted  to  execute  the  PSTOP3D  analysis.  A  shell  script  can  be  written 
to  interactively  generate  the  required  batch  script  as  shown  in  Appendix  A. 
A  typical  batch  script  to  execute  PSTOP3D  is  given  in  the  following  text. 
Green  text  is  text  the  user  must  tailor  to  their  particular  problem.  Red  text 
is  additional  explanation. 

! /bin/ksh 

#  pstop3d  qsub  script 

# 

#PBS  -N  ^ctp  C5| 

This  is  the  name  of  the  job.  Used  for  description  in  the  job  queue. 

#PBS  -1  select=B^: ncpus=32 :mpiprocs=32 

This  is  the  number  of  nodes  needed  on  Garnet.  Each  node  contains  32  processors. 
#PBS  -1  walltimeHOS  :  30  :  00| 

This  is  the  execution  time  in  hours : minutes : seconds . 

#PBS  -q  [standard  sMj 

This  is  the  name  of  the  standard  queue  for  small  problems. 

#PBS  -A  IkDCACCQUI^ 

This  is  the  account  number  needed  that  provides  time  resources  for  running  the 
j  ob . 

#PBS  -m  be _ 

#PBS  -M  [j  oe  .  somebody@usace  .  army .  milj 
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This  is  your  email.  An  email  will  be  sent  when  the  job  starts  running, 
cd  /u/|lserid/data/mctp  C5| 

Change  to  the  directory  where  the  PST0P3D  files  are  located, 
mkdir  /work/|aserid/mctp  ^ 

Make  a  work  directory. 

cp  /u/|lserid/data/mctp  C5|/^  /work/jlserid/mctp 

Copy  all  files  from  the  user  directory  to  the  work  directory, 
cd  /work/|nserid/mctp  C^ 

Change  to  the  work  directory.  It  is  important  to  change  the  directory  before  the 
remaining  lines  are  executed. 

cp  COMMAND.  FIL.fc^tc^.  txt  COMMAND.  FIL 

Copy  the  COMMAND. FIL  for  the  particular  problem  to  the  general  name  COMMAND. FIL. 
cp  COMMAND.  AUX.Inctp  aux  l| .  legacy .  txt  COMMAND.  AUX 

Copy  the  COMMAND. AUX  file  for  the  particular  problem  to  the  general  name 
COMMAND . AUX . 

cp  [ctJ_geo  .  dat .  0  .  origtopo  [ctJ__geo  .  dat .  0 . 2Douttopo 

Copy  the  topography  file  used  to  define  the  topography  for  2D  slice  output.  The 
.origtopo  file  is  the  original  unmodified  topography.  The  .2Douttopo  file  is  the 
modified  topography  that  includes  the  tops  of  the  buildings.  This  line  is  optional 
depending  on  whether  the  output  is  desired  at  a  distance  above  the  topography  or 
above  the  topography  and  tops  of  the  buildings. 

aprun  -n  |512|  .  /fe|stop3d  mctp  C5| 

This  line  submits  the  batch  file.  512  is  the  total  number  of  processors  (number  of 
nodes  x  32) .  pstop3d_mctp_C5  is  the  name  of  the  executable. 

exit 

To  submit  the  batch  script  to  the  job  queue,  type: 
qsub  batch_script_name 

where  “batch_script_name”  is  the  name  of  the  batch  script. 


ERDC  TR-15-5 


125 


8  Post-process  PSTOP3D  Results 

8.1  General 

The  information  contained  in  this  chapter  corresponds  to  Step  8  on  the 
process  diagram  shown  in  Figure  1.4.  Post-processing  of  the  results  from 
PSTOP3D  can  be  performed  using  a  MATLAB  utility  called  SNAPPS. 
SNAPPS  requires  the  use  of  several  files.  Three  of  the  files  are  produced  by 
PSTOP3D  and  will  need  to  be  downloaded  from  the  HPC  machine  to  the 
local  machine  where  MATLAB  is  installed.  The  files  needed  are: 

1.  RUNID_PX.SNP 

2.  RUNID_PY.SNP 

3.  RUNID_PZ.SNP 

4.  RUID.INFO 

The  RUNID  part  of  the  file  names  listed  above  is  explained  in  Section  7.6. 
The  other  files  required  pertain  to  the  topography  and  the  load  definition. 
These  will  be  requested  by  SNAPPS. 

The  .SNP  files  or  SNAP  files  are  the  output  files  produced  by  PSTOP3D 
that  have  been  temporally  decimated.  That  is,  the  files  were  produced 
using  the  INTERVAL_SNAP  discussed  in  Section  7.8.  The  _PX  file 
contains  the  results  for  a  section  cut  along  the  x  axis,  the  _PY  file  contains 
the  results  for  a  section  cut  along  the  y  axis,  and  the  _PZ  file  contains  the 
results  for  a  section  cut  along  the  z  axis.  The  location  of  the  section  is 
determined  by  the  IX_SNAP_PLN,  IY_SNAP_PLN,  and 
IZ_SNAP_PLN  data  items  explained  in  Section  7.8. 

8.2  How  to  execute  SNAPPS 

Once  the  output  files  have  been  downloaded  from  the  HPC  machine, 
SNAPPS  can  be  executed.  To  execute  SNAPPS: 

1.  Execute  MATLAB  in  Windows. 

2.  On  the  main  MATLAB  screen  shown  in  Figure  6.1,  navigate  to  the 
directory  where  the  scripts  are  installed.  The  main  script  to  execute  is 
called  snappsi53.m.  The  numbers  in  the  file  name  indicate  the  version 
number. 
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3.  The  path  may  need  to  be  changed  to  include  the  directory  and  sub¬ 
directories  where  the  scripts  are  installed.  On  the  Home  tab,  click  the  Set 
Path  button. 

4.  As  shown  in  Figure  6.2,  click  Add  with  Subfolders  and  select  the 
appropriate  directory.  Click  Save. 

5.  Double-click  on  the  snappsiss.m  script. 

6.  The  editor  with  snappsi53.m  will  appear.  Click  on  the  Run  button  on 
the  main  toolbar. 

7.  The  main  menu  for  the  processing  routines  will  appear  as  shown  in 
Figure  8.1. 

8.  On  the  menu,  select  the  button  View  or  extract  HPC  data. 

9.  The  dialog  in  Figure  8.2  will  appear.  Click  on  Specify  an  HPC  case 
("INFO")  file.  Navigate  to  where  the  INFO  file  is  saved  and  click  on  the 
file  name. 

10.  The  dialog  in  Figure  8.3  will  appear.  Click  Yes  to  allow  SNAPPS  to  change 
the  directory  to  where  the  data  file  is  located. 

11.  A  file  dialog  will  appear  for  the  input  of  the  MODEL_geo.dat.4  file  for 
these  data.  The  MODEL  parameter  is  explained  in  Section  7.6.  Navigate 
to  the  location  of  this  file  and  select  it.  This  file  is  located  in  a  directory 
called  legacy. 

12.  A  file  dialog  will  appear  for  the  input  of  the  RUN_ls.dat  file  that  was 
generated  by  the  rvloadset_io.m  script  as  explained  in  Section  6.6. 

13.  The  dialog  in  Figure  8.4  will  appear  to  select  the  section  cut  for  processing 
results.  This  corresponds  to  the  sections  defined  in  the  COMMAND.FIL 
file  described  in  Section  7.8. 

14.  The  dialog  in  Figure  8.5  will  appear  to  select  the  type  of  results  to  view. 
This  menu  applies  to  each  section  cut  selected  in  Step  13. 
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Figure  8.1.  Main  menu  for  SNAPPY. 


Figure  8.2.  Choose  data  to  visuaiize. 


Figure  8.3.  Specify  if 
SNAPPS  shouid  change 
or  create  directories. 

MENul  ^ 

Allow  Snapps  to 
change  and/or  create 
directories? 


Yes 


No 


Info 
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Figure  8.4.  Choose  section 
cut  to  extract  data. 

g  MENU  I  I  i~k 


Scope  of  Output  (plot  or  data 
extraction} 


Acoustic  in  xy -plane 


Figure  8.5.  Select  type  of  results  to 
view. 


8.3  Produce  animations  of  pressures 

To  produce  an  animation  of  the  pressures  over  time: 

1.  Select  either  the  All  outputs  (animation/single  frames)  or  Regions 
(animation/single  frames)  option  from  the  menu  shown  in  Figure  8.5. 
The  difference  in  the  options  is  the  All  outputs  will  output  for  the  entire 
analysis  region  and  the  Regions  will  allow  the  user  to  window  in  on  a 
specific  area. 

2.  If  Regions  are  selected,  the  dialog  in  Figure  8.6  will  request  the  region  to 
display.  Input  the  lower  left  and  upper  right  coordinates  of  the  box 
defining  the  region  to  plot. 
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3.  The  time  extents  will  be  requested  as  shown  in  Figure  8.7.  The  time 
extents  of  the  animation  can  be  specified  in  frames,  seconds,  or 
milliseconds.  For  each  option,  a  starting  and  ending  point  are  specified. 

For  the  data  shown  in  Figure  8.7,  First  frame  is  1  and  Last  frame  is 
410.  Since  there  are  410  output  frames,  the  animation  will  encompass  the 
entire  analysis  time  period. 

4.  The  color  limits  are  selected  as  shown  in  Figure  8.8.  The  symmetrical  color 
limit  will  cause  the  color  to  be  centered  around  o  Pa,  and  the  upper  and 
lower  limits  to  be  1  and  -1,  respectively. 

5.  The  dialog  shown  in  Figure  8.9  will  display  requesting  the  number  of 
contours  to  use  for  the  topography.  If  the  topography  has  been  modified  as 
discussed  in  Section  5.7,  then  the  contour  lines  will  include  the  tops  of  the 
buildings  as  shown  in  Figure  8.10.  Click  OK  to  re-plot  with  the  specified 
contours.  Click  Done  to  continue. 

6.  The  animation  will  be  constructed  and  displayed.  An  example  of  one  frame 
of  an  animation  is  shown  in  Figure  8.11. 

7.  After  the  animation  is  completed,  the  dialog  in  Figure  8.12  is  displayed. 
Options  are  available  to  regenerate  the  animation,  save  the  animation, 
modify  the  colors  or  contour  levels,  or  generate  a  Google  Earth  placemark. 

Figure  8.6.  Input  of  region. 
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Figure  8.7.  Selection  of  frames  for  animation. 


Figure  8.8.  Specify  color  limits  for  the  animation. 


I  Animation  parameters 


Specify  animation  details. 


The  colorscale  in  the  animation  can  be  determined  automatically  as  a 
percentage  of  the  maximum  output  amplitude  during  the  animation  period 
or  it  can  be  fixed  with  user-defined  limits. 


B  Automatic  color  limits: 

Saturate  colorscale  at  pet  of  max  output:  _ 1 


O  Specified  fixed  color  limits  (Pa) 
Max  colorscale  limit: 

Min  colorscale  limit: 


Symmetrical  color  limits  centered  on  0  (Pa): 

Colorscaie  limit:  1 


1 


-1 


Cancel 


OK 
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Figure  8.9.  Selection  of  number  of  contours. 


Figure  8.10.  Contour  of  topography. 
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Figure  8.11.  Example  of  frame  from  animation  of  pressures. 
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Figure  8.12.  Options  after  creation  of 
animation. 


8.4  Output  time  series  information 

There  were  several  locations  where  instrumentation  had  been  installed  on 

the  SMU  Campus.  At  these  locations,  the  time  histories  of  the  pressure 

variations  were  extracted.  The  locations  for  this  study  are  shown  in 

Figure  8.13.  To  plot  the  time  histories  at  specific  locations: 

1.  Convert  all  UTM  coordinates  into  local  coordinates.  To  perform  this 
conversion,  the  lower  left-hand  UTM  coordinates  will  be  subtracted  from 
each  of  the  desired  locations.  This  will  result  in  a  local  coordinate  system 
with  x=o,  and  y=o  in  the  lower  left-hand  comer. 

2.  From  the  menu  shown  in  Figure  8.5,  select  either  Time  series  at  grid 
points  or  Time  series  at  specific  points. 

3.  If  the  lime  series  at  grid  points  is  selected,  the  dialog  in  Figure  8.14  is 
displayed.  Input  the  lower  left  (xi,  yi)  and  upper  right  (x2,  y2)  comer  of 
the  region  to  extract  results.  Click  OK. 

4.  The  dialog  in  Figure  8.15  will  display.  Input  values  for  the  Number  of  x 
points  in  the  grid  and  Number  of  y  points  in  the  grid.  This  will 
form  a  grid  for  the  extraction  of  values. 

5.  Select  Centered  in  equal  area  regions  for  the  Locate  the  output 
points  option.  This  will  select  the  location  that  is  at  the  center  of  the  grid 


ERDC  TR-15-5 


133 


defined.  This  option  may  also  be  set  to  extract  results  at  the  nodes  of  the 
grid  defined. 

6.  The  time  extents  will  be  requested  as  shown  in  Figure  8.16.  The  time 
extents  of  the  time  history  can  be  specified  in  frames,  seconds,  or 
milliseconds.  For  each  option,  a  starting  and  ending  point  are  specified. 
For  the  data  shown  in  Figure  8.16,  First  frame  is  1  and  Last  frame  is 
410.  Since  there  are  410  output  frames,  the  time  history  will  encompass 
the  entire  analysis  time  period. 

7.  The  filter  dialog  as  shown  in  Figure  8.17  will  display.  For  the  version  of 
SNAPPS  used,  the  filter  functionality  was  not  fully  developed.  Therefore, 
check  the  No  filter  option  and  click  OK. 

8.  The  time  histories  of  the  pressures  are  output  at  the  selected  grid  locations 
as  shown  in  Figure  8.18. 

9.  An  option  will  be  displayed  to  save  these  histories  with  a  KML  if  desired. 

10.  The  time  histories  can  be  saved  as  a  MATLAB  time  series  object  as  shown 
in  Figure  8.19,  which  can  be  later  analyzed  with  the  Time  Series  Tools  in 
MATLAB. 


Figure  8.13.  Locations  for  output  of  time  history  information. 
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Figure  8.14.  Select  region  to  compute  time  series  information. 


Figure  8.15.  Select  grid  size  for  time  series  information. 


ERDC  TR-15-5 


135 


Figure  8.16.  Select  the  extents  of  the  time  history. 


Figure  8.17.  Select  filter  for  results. 


The  Nyquist  frequency  for  this  data  set  is  20  Hz.  The  data  set  will  be  bandpassed 
according  to  the  following  specification.  Frequencies  specified  in  excess  of  the 
Nyquist  will  be  set  1  percent  below  the  Nyquist.  Minimum  allowable  frequency  is 
0.1  Hz.  Stop  and  pass  limits  cannot  be  specified  as  equal. 


O  Browse  to  an  existing  filtered  data  set 


O  New  filter 

Bandpass  lower  stop  limit  (Hz):  _ 0.1 

Bandpass  lower  pass  limit  (Hz):  _ 0.3 

Bandpass  upper  pass  limit  (Hz):  _ 19.6 

Bandpass  upper  stop  limit  (Hz):  _ 19.8 

CutdB:  I  60~ 

Ripple  dB:  |  1 

Filter  type:  Butterworth 

O  View  the  filter  design 


No  filter 


Cancel 


OK 
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Figure  8.18.  Time  history  output  at  grid  points. 


2468  10  2468  10 

t(s) 


Figure  8.19.  Prompt  to 
save  time  history  signais. 


To  enter  the  locations  at  specific  points: 

1.  From  the  menu  shown  in  Figure  8.5,  select  Time  series  at  specific 
points. 

2.  The  dialog  in  Figure  8.20  will  display.  Select  the  Point-by-point  option. 

3.  The  dialog  in  Figure  8.21  will  display.  Input  the  x  and  y  location  of  a  point 
and  click  the  Add  this  location  button. 

4.  Continue  to  add  points  and  click  the  Add  this  location  button  each  time. 

5.  For  the  last  point,  click  Add  this  location  and  then  click  Done. 

6.  Select  the  time  extents  as  shown  in  Figure  8.16  and  the  filter  option  as 
shown  in  Figure  8.17. 
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7.  The  dialog  in  Figure  8.22  will  display.  Select  Plot  the  signal  using 
format  for  individual  time  history.  Each  input  location  will  be 
plotted  as  shown  in  Figure  8.23. 

8.  If  the  option  Plot  the  signal  using  format  for  multi-signal  on 
common  axis  is  selected,  the  time  histories  will  be  plotted  as  shown  in 
Figure  8.24. 

9.  Click  Save  this  signal  to  a  file  in  Figure  8.22  to  save  the  time  history  to 
a  file.  The  options  shown  in  Figure  8.19  will  display. 

Figure  8.20.  Specify  output  at 
points. 


Figure  8.21.  input  of  points  at  which  to  output  time 
histories. 
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Figure  8.22.  Select  how  to  plot  signals. 


Figure  8.23.  Time  history  of  pressure  at  single  point. 


t(t» 


Figure  8.24.  Time  histories  of  pressure  at  multiple  points. 
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3635390.344 
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Boaz  -  East,  Channel  2  -  L-4 
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1129.376 
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Collins,  Channel  2 -L-4 

4 
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3 

707431.968 
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1 
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If  the  time  histories  are  saved  to  a  MATLAB  time  series  object  file,  then 

the  MATLAB  time  series  tools  can  be  used  to  display  and  perform  data 

manipulations.  To  read  and  display  a  time  series  object  file: 

1.  Execute  MATLAB. 

2.  Navigate  to  the  directory  where  the  file  is  saved  using  the  Current 
Folder  window  in  MATLAB.  The  file  will  have  a  .mat  extension. 

3.  Right-click  on  the  file  and  select  Load. 

4.  The  time  series  object  will  appear  in  the  Workspace  window  of 
MATLAB. 

5.  Right-click  the  time  series  object  and  select  Open  in  Time  Series  Tools. 

6.  Select  the  type  of  plot,  e.g..  Spectral  Plots  or  Time  Plots,  as  shown  in 
Figure  8.25.  Click  Display. 

7.  Both  the  spectral  and  time  plots  are  shown  in  Figure  8.26.  The  spectral 
plots  indicate  that  for  the  time  histories  plotted,  the  dominant  frequencies 
are  1.8  and  3  Hz. 


Figure  8.25.  Time  series  toois  in  MATLAB. 
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Figure  8.26.  Time  and  spectral  plots  in  MATLAB. 


8.5  Plot  peak  pressures 

The  peak  pressures  can  be  plotted  for  the  entire  analysis  region  as  follows: 

1.  Select  Peak  values  map  from  the  menu  shown  in  Figure  8.5. 

2.  The  time  extents  will  be  requested  as  shown  in  Figure  8.27.  The  time 
extents  of  the  peak  values  can  be  specified  in  frames,  seconds,  or 
milliseconds.  For  each  option,  a  starting  and  ending  point  are  specified. 
For  the  data  shown  in  Figure  8.27,  First  frame  is  1  and  Last  frame  is 
410.  Since  there  are  410  output  frames,  the  peak  values  will  encompass 
the  entire  analysis  time  period. 

3.  The  peak  values  for  the  x-y  plane  are  plotted  as  shown  in  Figure  8.28. 
Roads  and  buildings  have  been  added  to  the  plot  with  Global  Mapper. 
Refer  to  Section  8.8  for  instructions  on  how  to  add  GIS  data  to  a  results 
plot.  An  example  of  peak  pressures  for  the  y-z  plane  is  shown  in 
Figure  8.29. 

4.  The  dialog  shown  in  Figure  8.9  will  display  requesting  the  number  of 
contours  to  use  for  the  topography.  If  the  topography  has  been  modified  as 
discussed  in  Section  5.7,  then  the  contour  lines  will  include  the  tops  of  the 
buildings.  The  Number  of  contour  levels  and  the  Color  are  input. 
Click  OK  to  re-plot  with  the  specified  contours.  Click  Done  to  continue. 
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5.  After  the  peak  values  are  plotted,  the  dialog  in  Figure  8.30  is  displayed. 
Options  are  available  to  save  the  image  as  a  KML  placemark,  modify  the 
contours  and  colors,  or  compare  to  a  reference  set  of  pressures.  The 
Compare  to  reference  option  is  explained  in  Section  8.7. 

Figure  8.27.  Specify  time  extents  for  peak 
vaiues. 


Figure  8.28.  Peak  pressures  in  x-y  piane  for  5  Hz  source  in  upper  right 
corner  for  410  frames. 
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Figure  8.29.  Peak  pressures  in  y-z  piane  for  source  iocated  on  interior  for 

410  frames. 


Figure  8.30.  Options  for 
peak  vaiues. 


The  peak  pressure  plot  gives  the  peak  pressure  at  every  nodal  point 
location  in  the  analysis  region  over  the  time  period  specified.  If  the  first 
frame  value  in  Figure  8.27  is  input  as  200  and  the  last  frame  value  is  200, 
the  peak  pressures  over  the  entire  region  are  computed  at  each  nodal  point 
for  frame  200  (or  5  sec)  as  shown  in  Figure  8.31.  If  the  first  frame  value  in 
Figure  8.27  is  input  as  1  and  the  last  frame  value  is  200,  the  peak  pressures 
over  the  entire  region  are  computed  at  each  nodal  point  for  frames  1 
through  200  (or  o  sec  through  5  sec)  as  shown  in  Figure  8.32. 
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Figure  8.31.  Peak  pressures  for  frame  200  (5  sec) 
for  source  in  upper  ieft  corner. 


Figure  8.32.  Peak  pressures  for  frame  1  through 
200  (5  sec)  for  source  in  upper  ieft  corner. 


8.6  Produce  integrated  energy  plots 

The  integrated  energy  can  be  plotted  for  the  entire  analysis  region  as 

follows: 

1.  Select  Integrated  energy  map  from  the  menu  shown  in  Figure  8.5. 

2.  The  time  extents  will  be  requested  as  shown  in  Figure  8.33.  The  time 
extents  of  the  integrated  energy  can  be  specified  in  frames,  seconds,  or 
milliseconds.  For  each  option,  a  starting  and  ending  point  are  specified. 
For  the  data  shown  in  Figure  8.33,  First  frame  is  1  and  Last  frame  is 
410.  Since  there  are  410  output  frames,  the  integrated  energy  will 
encompass  the  entire  analysis  time  period. 


ERDC  TR-15-5 


144 


3.  The  integrated  energy  for  the  x-y  plane  is  plotted  in  Figure  8.34.  The  red 
lines  in  the  figure  are  contours  of  the  topography.  The  solid  red  areas  are 
buildings. 

4.  The  dialog  shown  in  Figure  8.9  will  display  requesting  the  number  of 
contours  to  use  for  the  topography.  If  the  topography  has  been  modified  as 
discussed  in  Section  5.7,  then  the  contour  lines  will  include  the  tops  of  the 
buildings.  The  Number  of  contour  levels  and  the  Color  are  input. 
Click  OK  to  re-plot  with  the  specified  contours.  Click  Done  to  continue. 

5.  After  the  integrated  energy  values  are  plotted,  the  dialog  in  Figure  8.30  is 
displayed.  Options  are  available  to  save  the  image  as  a  KML  placemark, 
modify  the  contours  and  colors,  or  compare  to  a  reference  set  of  pressures. 
The  Compare  to  reference  option  is  explained  in  Section  8.7. 

6.  The  integrated  energy  plot  gives  the  integrated  energy  at  every  nodal  point 
location  in  the  analysis  region  over  the  time  period  specified.  If  the  first 
frame  value  in  Figure  8.27  is  input  as  200  and  the  last  frame  value  is  200, 
the  integrated  energy  over  the  entire  region  is  computed  at  each  nodal 
point  for  frame  200  (or  5  sec)  as  shown  in  Figure  8.35.  If  the  first  frame 
value  in  Figure  8.27  is  input  as  1  and  the  last  frame  value  is  200,  the 
integrated  energy  over  the  entire  region  is  computed  at  each  nodal  point 
for  frames  1  through  200  (or  o  sec  through  5  sec)  as  shown  in  Figure  8.36. 

Figure  8.33.  Integration  period  for  energy. 
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Figure  8.34.  Integrated  energy  for  the  x-y  plane  for  a  source  in  the  upper 

right  corner. 


Figure  8.35.  Integrated  energy  for  frame  200  (5  ec)  for  source  in 
upper  left  corner. 
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Figure  8.36.  Integrated  energy  for  frame  1  through  200  (5  sec)  for 
source  in  upper  left  corner. 


8.7  Compare  to  reference 

The  Compare  to  reference  option  available  for  peak  pressures  and 
integrated  energy  allows  the  comparison  of  a  set  of  values  to  a  reference 
set.  The  comparison  is  made  by  subtracting  the  reference  set  from  the 
current  values  plotted.  This  allows  the  difference  between  the  two  sets  of 
values  to  be  highlighted.  To  specify  a  reference  set: 

1.  After  a  set  of  values  have  been  plotted,  the  menu  shown  in  Figure  8.37  is 
displayed.  The  menu  lists  several  options  for  selecting  the  reference  set. 

2.  As  the  plots  are  made  in  MATLAB,  they  are  also  saved  to  disk  with  a  .fig 
extension.  If  the  Browse  to  .fig  file  option  is  selected,  a  previous  plot 
maybe  selected  as  the  reference  source. 

3.  Browse  to  the  stored  figure  and  select  the  file.  The  dialog  shown  in  Figure 
8.38  is  displayed.  Select  the  number  of  colors  and  the  min  and  max  values. 
Click  OK  The  comparison  plot  will  be  generated. 

Figure  8.39  shows  the  peak  pressures  for  a  10  Hz  source  in  the  upper  right 
corner  of  the  urban  area.  If  Figure  8.28,  which  shows  the  peak  pressures  for 
a  5  Hz  source,  is  used  as  the  reference,  the  results  are  shown  in  Figure  8.40. 
The  results  indicate  that  the  sunken  highway  is  acting  as  a  waveguide.  The 
purple  color  indicates  an  increase  from  the  reference  values. 
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Another  example  of  comparing  to  a  set  of  reference  values  is  shown  in 
Figure  8.41  through  8.43.  In  Figure  8.41,  the  peak  pressures  for  a  5  Hz 
source  in  the  upper  right  corner  of  the  analysis  region  are  shown. 

Figure  8.42  shows  the  peak  pressures  for  the  same  source,  but  with  the 
buildings  aggregated  into  larger  areas.  The  difference  plot  shown  in 
Figure  8.43  shows  that  in  the  region  of  the  sunken  highway,  the  aggregation 
did  not  make  a  significant  difference. 


Figure  8.37.  Compare 
a  set  of  values  to  a 
reference  set. 


Figure  8.38.  Options  for  reference  plot. 
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Figure  8.39.  Peak  pressures  in  x-y  plane  for  10  Hz  source  in  upper  right  corner  for 

410  frames. 
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Figure  8.40.  Results  of  reference  calculation  using  5  and  10  Hz 

sources. 
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Figure  8.41.  Peak  pressures  for  5  Hz  source  located  in  upper 
right  corner  for  individual  buildings. 


Figure  8.42.  Peak  pressures  for  5  Hz  source  located  in  upper 
right  corner  for  aggregated  buildings. 
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Figure  8.43.  Comparison  plot  for  difference  between  individual 
buildings  and  aggregated  buildings. 


8.8  Using  Global  Mapper  to  enhance  visual  display  of  results 

Graphics  of  the  results  that  are  generated  by  the  SNAPPS  utility  in 

MATLAB  can  be  enhanced  for  presentation  by  combining  the  results  with 

GIS  features.  This  can  be  performed  using  the  Global  Mapper  program. 

For  example,  the  difference  plot  shown  in  Figure  8.40  was  enhanced  by 

adding  roads  and  buildings  to  the  plot.  To  enhance  a  picture: 

1.  Execute  Global  Mapper  and  have  the  GIS  information  for  the  area 
displayed. 

2.  Navigate  to  the  directory  where  the  plot  is  saved  that  is  to  be  enhanced. 
The  plot  can  be  in  several  bitmap  formats  including  BMP,  JPG,  and  TIFF. 

3.  Click  on  the  file  name  in  Windows  Explorer  and  drag  the  file  to  the  Global 
Mapper  window. 

4.  Global  Mapper  will  ask  if  the  user  would  like  to  rectify  the  image.  Select 
Yes. 

5.  Global  Mapper  loads  the  image  and  displays  it  as  shown  in  Figure  8.44.  In 
the  leftmost  window,  click  and  drag  a  box  to  zoom  to  an  area  to  select  the 
first  point  to  rectify. 

6.  In  the  middle  window,  the  zoomed-in  area  is  displayed.  Click  the  mouse  to 
select  the  upper-left  point  of  the  region. 
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7.  Zoom  to  the  same  area  in  the  right-most  window.  Click  to  select  the  same 
location  of  the  upper  left  corner  on  the  GIS  data.  Click  the  Add  Point  to 
List  button. 

8.  Continue  to  zoom  and  select  the  comers  of  the  box  defining  the  analysis 
region.  Only  four  points  will  be  required  to  rectify  the  image. 

9.  After  the  final  and  fourth  point  is  selected  as  shown  in  Figure  8.45,  click 
the  OK  button. 

10.  The  image  will  be  rectified  to  the  GIS  area  and  displayed.  Now  any  GIS 
information  can  be  overlaid  upon  the  image. 

11.  The  completed  figure  with  roads  and  buildings  added  is  shown  in 
Figure  8.46. 

Using  Global  Mapper,  the  rectified  image  can  have  any  of  the  GIS  data 
overlaid  upon  the  image  to  enhance  presentation.  The  image  can  be  saved 
from  Global  Mapper  in  a  variety  of  image  formats  for  insertion  in 
Microsoft  Word  or  PowerPoint.  The  image  can  also  be  copied  from 
desktop  by  using  the  CTRL- ALT  key  sequence  and  pasting  directly  into 
Word  or  PowerPoint  where  it  can  be  cropped  to  the  desired  size. 


Figure  8.44.  Selecting  point  1  to  rectify  an  image  to  GIS  data. 
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Figure  8.45.  Selecting  point  4  to  rectify  an  image  to  GIS  data. 


Figure  8.46.  Final  rectified  figure  with  buildings  and  roads  added. 
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9  Conclusions  and  Implications 

The  PSTOP3D  finite  difference  time  domain  program  can  be  used  to 
perform  an  acoustical  analysis  of  a  3-D  urban  environment  in  the 
infrasound  passband.  The  modeling  for  the  environment  can  be 
accomplished  using  the  commercial  GIS  programs  Global  Mapper  and 
ArcMAP.  Items  such  as  roads,  streets,  highways,  parking  lots,  and 
buildings  can  be  modeled  to  effectively  represent  the  urban  environment. 

The  results  from  the  3-D  urban  models  allow  users  to  understand 
propagation  path  effects  in  an  urban  setting.  The  output  of  the  model  is  a 
peak  energy  map  of  the  area,  which  could  be  used  in  pre-deployment 
planning  for  infrasound  array  installation  selection  to  determine  the 
propagation  effects  expected  at  a  given  location.  Results  of  a  real-world 
modeled  environment  compared  to  infrasound  and  seismic  data  collected 
will  be  discussed  in  a  future  report  in  this  report  series. 
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Appendix  A 

This  appendix  contains  a  Unix  script  that  generates  a  batch  script  for 
running  a  PSTOP3D  analysis.  This  script  simplifies  the  running  of 
PSTOP3D  by  taking  input  from  the  user  and  generating  a  script  for  the 
specific  problem  being  analyzed.  The  script  is  as  follows: 

# ! /bin/ksh 

#  Write  runpstop.sh 

#  754=rwxr-xr--  Group  can  read  and  execute 

#  Use  "chmod  754  pstop3d_writeQsubScript_ksh . sh"  if  not  executable 

#  To  use,  current  directory  must  have  this  script  and  all  pstop3d  input  files 

#  Current  directory  will  be  the  output  directory 

#  - 

#  Print  number  of  arguments 
echo  "Number  of  arguments  is  $#" 

#  Test  for  correct  number  of  arguments 
if  (test  "$#"  =  "1")  then 

echo  >>  /dev/null 
elif  (test  "$#"  =  "4")  then 
echo  >>  /dev/null 
elif  (test  "$#"  =  "6")  then 
echo  >>  /dev/null 
elif  (test  "$#"  =  "7")  then 
echo  >>  /dev/null 
else 
echo 

echo  "The  correct  usage  is:" 

echo  "  ./makerun.sh  datafile  [nodes  walltime  queue]  [topo  auxfile]  [account]" 
echo  "  datafile  =  data  file  name  without  extension." 

echo  "  Identifier  for  the  set  of  data  files,  e.g.,  COMMAND.FIL.identifier.txt" 
echo  "  nodes  =  number  of  nodes  to  use  (default  is  4) ." 

echo  "  walltime  =  time  limit  of  run  in  HH:MM:SS  (default  is  00:30:00)  ." 

echo  "  queque  =  debug  or  standard  (default  is  standard) . " 

echo  "  topo  =  topographic  surface  to  use  for  slice  output." 

echo  "  orig  =  for  original  topography." 

echo  "  mod  =  for  modified  topography  (default) " 

echo  "  auxfile  =  auxiliary  file  to  use  from  COMMAND.AUX.identifier.legacy.txt" 
echo  "  account  =  account  to  run  job  under  (default  is  ERDCNUM) ." 
echo 
exit 
fi 

#  Set  defaults 
dfile=$l 
nodes=4 

timelimit=" 00 : 30 : 00" 
qname=" standard_sm" 
topo="mod" 
account="ERDCNUM" 
auxf ile=$df ile 

#assign  values  from  parameters  entered 

#  Test  arguments  for  queue  name, 
if  (test  "$#"  =  "4")  then 

if  (test  "$4"  =  "debug")  then 
echo  >>  /dev/null 
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elif  (test  "$4"  =  " standard_sm" )  then 

echo  >>  /dev/null 

else 

echo 

echo  "Invalid  queue  name  entered.  Valid  names  are  debug  or  standard_sm. " 
exit 
fi 
fi 

if  (test  "$#"  =  "4")  then 
nodes=$2 
timelimit=$3 
qname=$4 
fi 

if  (test  "$#"  =  "6")  then 
nodes=$2 
timelimit=$3 
qname=$4 
topo=$5 
auxfile=$6 
fi 

if  (test  "$#"  =  "7")  then 
nodes=$2 
timelimit=$3 
qname=$4 
topo=$5 
auxf ile=$  6 
account=$7 
fi 

#  Define  directories 

#  Customize  for  your  individual  directory. 
thisdir= 'pwd' 

j  obid=$ { df ile }_$$ 
workdir=$WORKDIR/ $thisdir 

#  Create  run  file 

runFile="run_pstop3d_$ { df ile }_by_qsub_$nodes . txt" 

echo  "#!/bin/ksh"  >  $runFile 

echo  "#  pstop3d  qsub  script"  >>  $runFile 

echo  "#"  >>  $runFile 

#  Command  line  OK.  Start  job. 

#  Add  PBS  lines  using  required  arguments 
echo  "#PBS  -N  $dfile"  >>  $runFile 

echo  "#PBS  -1  select=$nodes : ncpus=32 : mpiprocs=32 "  >>  $runFile 
echo  "#PBS  -1  walltime=$timelimit"  >>  $runFile 
echo  "#PBS  -q  $qname"  >>  $runFile 
echo  "#PBS  -A  $account"  >>  $runFile 

#  Change  dir  to  current  directory  where  job  executed  from.  This  should  be  the 
directory 

#  that  contains  the  data  files.  Need  to  do  this  so  script  will  work  from  qsub. 
echo  "cd  $thisdir"  >>  $runFile 

echo  "mkdir  $workdir"  >>  $runFile 

#  Copy  needed  files 

echo  "cp  $thisdir/*  $workdir"  >>  $runFile 
echo  "cd  $workdir"  >>  $runFile 

#  Add  lines  to  copy  topography  to  use  the  .2Douttopo  file 
if  (test  "$topo"  =  "orig")  then 

echo  "cp  $ { df ile }_geo . dat . 0 . origtopo  $ { df ile }_geo . dat . 0 . 2Douttopo"  »  $runFile 
fi 

#  Add  lines  to  copy  files  to  required  names 

#echo  'echo  " .  Copying  files  to  work  directory  -  $workdir" '  >>  $runFile 

echo  "cp  COMMAND. FIL. $df ile . txt  COMMAND. FIL"  »  $runFile 

echo  "cp  COMMAND. AUX.$auxf ile. legacy. txt  COMMAND. AUX"  »  $runFile 
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#  Add  lines  to  execute  code  and  exit 

#echo  'echo  " .  Submitting  batch  job"'  >>  $runFile 

echo  "aprun  -n  ( ( $nodes*32 ) )  . /pstop3d_$df ile"  >>  $runFile 

echo  "exit"  >>  $runFile 

#  Ask  to  run 

echo  "Do  you  want  to  run  now?  \c" 
read  answer 

if  (test  "$answer"  =  "y")  then 
# submit  job 
qsub  $runFile 
else 
exit 
fi 
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